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ABSTRACT 

Many parameters constraining the spectral appearance of exoplanets are still poorly understood. We therefore 
study the properties of irradiated exoplanet atmospheres over a wide parameter range including metallicity, C/0 
ratio and host spectral type. We calculate a grid of 1-d radiative-convective atmospheres and emission spectra. 
We perform the calculations with our new Pressure-Temperature Iterator and Spectral Emission Calculator for 
Planetary Atmospheres (PETIT) code, assuming chemical equilibrium. The atmospheric structures and spectra 
are made available online. We find that atmospheres of planets with C/0 ratios ~ 1 and Teff > 1500 K can 
exhibit inversions due to heating by the alkalis because the main coolants CH 4 , H 2 O and HCN are depleted. 
Therefore, temperature inversions possibly occur without the presence of additional absorbers like TiO and 
VO. At low temperatures we find that the pressure level of the photosphere strongly influences whether the 
atmospheric opacity is dominated by either water (for low C/0) or methane (for high C/0), or both (regardless 
of the C/0). Eor hot, carbon-rich objects this pressure level governs whether the atmosphere is dominated by 
methane or HCN. Eurther we find that host stars of late spectral type lead to planetary atmospheres which have 
shallower, more isothermal temperature profiles. In agreement with prior work we find that for planets with 
Teff < 1750 K the transition between water or methane dominated spectra occurs at C/0 ~ 0.7, instead of ~ 1, 
because condensation preferentially removes oxygen. 

Keywords: methods: numerical — planets and satellites: atmospheres — radiative transfer 


1. INTRODUCTION 

In a number of existing studies the range of possible C/0 
ratios in protoplanetary disks and the resulting implications 
for the C/0 ratios in the gaseous envelopes of extrasolar plan¬ 
ets is investigated (see, e.g., Oberg et al. 2011; Ali-Dib et al. 
2014; Helling et al. 2014; Marboeuf et al. 2014; Thiabaud 
et al. 2014). These kind of studies are interesting, as they 
may help to predict the spectral appearance of atmospheres 
of planets formed via different pathways in the circumstel- 
lar disks. In a further example, Madhusudhan et al. (2014a) 
studies the range of possible C/0 ratios for 2 different disk 
models, depending on the formation and migration mecha¬ 
nism invoked to form hot jupiters. The result of these studies 
is that large planetary C/0 ratios, close to unity, are possible 
even when considering disks of solar composition (the solar 
value is C/O© ~ 0.55, see Asplund et al. 2009). For disks 
with supersolar C/0 ratios the planetary C/0 ratios should be 
even higher, although stars with C/0 ratios close to and big¬ 
ger than 1 may be quite rare (Fortney 2012). The C/0 ratio 
is particularly interesting for the spectral appearance of exo¬ 
planets because for high enough temperatures (T > 1000 K) a 
C/0 < 1 giant planet will have appreciable amounts of H 2 O in 
its atmosphere and almost no CH 4 , whereas for C/0 > 1 the 
situation is the opposite and CH 4 is much more abundant than 
H 2 O. This transition happens quite sharply (see, e.g., Koppa- 
rapu et al. 2012; Madhusudhan 2012). Furthermore, conden¬ 
sation processes can potentially lead to local C/0 ratios of ~ 
1-2 in the gas phase, even if the global atmospheric C/0 ratio 
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is smaller than 1 (see Helling et al. 2014). The reason this is 
for the locking up of oxygen in silicates, as has already been 
suggested by Fortney et al. (2006). 

Both H 2 O and CH 4 have strong absorption features and 
their main absorption bands between -1.3 and 5 fim are alter¬ 
nately located in wavelength space. Thus hot gaseous planets 
with C/0 < 1 and C/0 > 1 in the spectrally active regions 
should be quite easily distinguishable and might give hints 
on the planet’s formation history such as the location of for¬ 
mation in the protoplanetary disk and its migration through 
it (Madhusudhan et al. 2014a). For even higher temperatures 
(T > 1750 K), and C/0 > 1, HCN takes over as the most im¬ 
portant carbon-carrying infrared absorber as it becomes more 
abundant than CH 4 in the spectrally active parts of the atmo¬ 
spheres (see, e.g., Moses et al. 2013). “Spectrally active” 
denotes the regions where the radiation seen in the planet’s 
emergent spectrum originates. The respective atmospheres 
are then not dominated by CH 4 anymore, but by HCN. Distin¬ 
guishing H 2 O and HCN absorption features should be possi¬ 
ble, due to the different spectral signatures of HCN and H 2 O 
in the NIR and IR. Therefore a distinction between O and C 
dominated atmospheres is possible also at high temperatures. 

Motivated by the fundamentally different spectral appear¬ 
ances of the two C/0 cases, Madhusudhan (2012) proposed 
a 2 -d classification scheme for characterizing giant extrasolar 
planets, using the C/0 ratio and the incident stellar flux as di¬ 
mensions. In his work, the importance of CH 4 for the C/0 > 
1 cases is most strongly emphasized, but the possible impor¬ 
tance of HCN and C 2 H 2 is mentioned as well. A 1-d classi¬ 
fication scheme for hot giant planets, based only on the stel¬ 
lar flux, had already been proposed by Fortney et al. (2008) 
before, featuring “cold” hot jupiters without a temperature in¬ 
version and “hot” hot jupiters with a temperature inversion 
caused by TiO and VO absorption. However, some “hot” hot 
jupiters are not thought to have a inversion, contradicting the 
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Opacity source 

Spectral range [jum] 

Line list 

Partition function 

Pressure broadening 

CH 4 

0.83 < A 

Yurchenko and Tennyson (2014) 

(a) 

(a) 

CH 4 

0.86 < A 

Rothman et al. (2013) 

Fischer et al. (2003) 

7 air. Rothman et al. (2013) 

C 2 H 2 

1<A< 16.5 

Rothman et al. (2013) 

Fischer et al. (2003) 

7 air, Rothman et al. (2013) 

CO 

1.18<T 

Rothman et al. (2010) 

Fischer et al. (2003) 

7 air, Rothman et al. (2010) 

CO 

0.112<T<0.43 

Kurucz (1993) 

Fischer et al. (2003) 

Eq. (15), Sharp and Burrows (2007) 

CO 2 

1 < T < 38.76 

Rothman et al. (2010) 

Fischer et al. (2003) 

7 air, Rothman et al. (2010) 

H 2 S 

0.88 < A 

Rothman et al. (2013) 

Fischer et al. (2003) 

7 air, Rothman et al. (2013) 

H 2 

0.28 < A 

Rothman et al. (2013) 

Fischer et al. (2003) 

7 air, Rothman et al. (2013) 

H 2 

0.08<T<0.18 

Kurucz (1993) 

Fischer et al. (2003) 

Eq. (15), Sharp and Burrows (2007) 

HCN 

2.92 < A 

Harris et al. (2006), 

Barber et al. (2014) 

Fischer et al. (2003) 

Eq. (15), Sharp and Burrows (2007) 

H 2 O 

0.33 < A 

Rothman et al. (2010) 

Fischer et al. (2003) 

7 air, Rothman et al. (2010) 

K 

0.05 < A 

Piskunov et al. (1995) 

Sauval and Tatum (1984) 

N. Allard, Schweitzer et al. (1996) 

Na 

0.1 <A 

Piskunov et al. (1995) 

Sauval and Tatum (1984) 

N. Allard, Schweitzer et al. (1996) 

NHs 

1.43 <T 

Rothman et al. (2013) 

Fischer et al. (2003) 

7 air, Rothman et al. (2013) 

OH 

0.52 < A 

Rothman et al. (2010) 

Fischer et al. (2003) 

7 air, Rothman et al. (2010) 

PH 3 

2.78 < A 

Rothman et al. (2013) 

Fischer et al. (2003) 

7 air, Rothman et al. (2013) 


Table 1 

References for the atomic and molecular opacities used in the PETIT code, (a): We use precalculated cross-sections from the ExoMol website. For these only 

Doppler broadening has been taken into account. 


1-d classification system. Madhusudhan (2012) argued that 
this could possibly be explained using the 2-d classification 
scheme, as TiO and VO should not be very abundant in plan¬ 
ets with a high C/0 ratio. In addition, there could be further 
reasons why TiO and VO should not be in the upper part of the 
atmosphere, e.g. due to settling and inefficient vertical mix¬ 
ing, cold-trap depletion or photodissociation (Spiegel et al. 
2009; Showman et al. 2009; Parmentier et al. 2013; Knutson 
et al. 2010). 

Observational evidence for planets with C/0 > 1 is scarce 
and the most prominent case, WASP-12b (Madhusudhan et al. 
2011), is controversial (Crossfield et al. 2012; Swain et al. 
2013; Stevenson et al. 2014). Current analyses of the photo¬ 
metric data indicate a C/0 ratio < 1: Line et al. (2014) esti¬ 
mated C/0 ratios for 9 hot jupiters and found that while in 7 
out of 9 cases (HD 209458b, GJ436b, HD 149026b, WASP- 
12b, WASP-19b, WASP-43b, TrES-2b) a C/0 value of 1 was 
within their 1 cr confidence interval, in 6 out of 9 cases the so¬ 
lar value was within the 1 cr interval as well. Benneke (2015) 
analyzed 8 hot jupiters (HD 209458b, WASP-19b, WASP- 
12b, HAT-P-lb, XO-lb, HD189733b, WASP-17b and WASP- 
43b) using a self-consistent retrieval analysis, ruling out C/0 
> 1 for all of them. Clearly the quality and quantity of the 
photometric and spectroscopic observations needs to improve 
before more conclusive results can be obtained for many of 
these planets (Line et al. 2013). 

A further example for a planet with a C/0 ratio close to 1 is 
HR 8799b, for which C/0 = 0.96 + 0.01 or 0.97+^-^^ has been 
estimated (Lee et al. 2013), depending on whether clouds are 
included in the model or not. 

Although all the C/0 ratios obtained by the above studies 
are depending on the assumptions made in the various re¬ 
trieval models, the current analysis of data does not indicate 
any planet with C/0 > 1. Further, while the current qual¬ 
ity of data is still too low for obtaining reliable retrieval re¬ 
sults in many cases, upcoming observing facilities such as 
JWST should greatly help to decipher the composition of hot 
jupiters. 

In conclusion, the C/0 ratio, together with the effective 
temperature, should be a key parameter constraining a hot 
jupiter’s spectral appearance and thus we want to study how 
the interplay between the C/0 ratio and other parameters af¬ 


fect the atmospheres. Systematic studies of exoplanet atmo¬ 
spheres have been published in the literature before (see, e.g., 
Sudarsky et al. 2003), and although the C/0 has been sug¬ 
gested to be of importance already a decade ago (Seager et al. 
2005), no systematic study of the atmospheric characteristics 
as a function of the C/0 ratio has been carried out so far. 
Therefore, we publish a grid of emission spectra and pressure 
temperature {PT) structures for self-consistent hot jupiters at¬ 
mospheres for varying C/0 ratio, [Fe/H], distance to the star, 
stellar host spectral type and planetary log(g). 

The results were calculated with our new Pressure- 
Temperature Iterator and Spectral Emission Calculator for 
Planetary Atmospheres (PETIT) code. PETIT solves the 1- 
d plane parallel structure of the atmosphere assuming local 
thermal equilibrium (LTE), radiative equilibrium or convec¬ 
tion and equilibrium chemistry. Our goal is to investigate 
the behavior of planetary atmospheres in the parameter space 
covered by our grid. Furthermore we make the atmospheric 
PT -structures, abundance profiles, and resulting spectra pub¬ 
licly available for use in, e.g., the evaluation of observational 
data of planetary emission spectra. 

In Section 2 we introduce our code and explain its individ¬ 
ual modules. We also show some of the tests we carried out to 
check the results of our code for consistency. In Section 3 we 
discuss how the assumptions in our code constrain the param¬ 
eter range of the atmospheric grid. In Section 4 we report on 
how the grid was set up and how the calculations were carried 
out. The results can be found in Section 5, the discussion and 
conclusions are in Section 6. 

2. DESCRIPTION OF THE CODE 
2.1. Opacity database 

The current version of our opacity database comprises 
atomic and molecular line and continuum opacities from 
ultra-violet to infrared wavelengths. So far only absorption 
processes are treated. Scattering of incoming stellar radia¬ 
tion by molecules in the planetary atmosphere causes the stel¬ 
lar photons to traverse, on average, a somewhat longer dis¬ 
tance through the atmosphere before reaching a certain pres¬ 
sure level. Hence, the photons will on average be absorbed 
at slightly lower pressures (higher altitudes) than if absorp¬ 
tion only is considered. Because the reported optical albedos 
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of hot jupiters are very low, in the low single digit percent¬ 
age range, as summarized by Madhusudhan et al. (2014b), 
absorption appears much more important than scattering in 
these objects. Therefore, include only absorption in the ra¬ 
diative transfer calculation One exception is mentioned, how¬ 
ever, with Kepler 7-b having a geometric albedo of 0.32 + 
0.03 (Demory et al. 2011). 

2.1.1. Molecular and atomic line opacities 

A list of all line opacity sources, together with a reference to 
the corresponding line lists, pressure broadening parameters 
and partition functions can be found in Table 1 . Our method 
to speed up molecular opacity calculations is explained in Ap¬ 
pendix A. 

All molecular and atomic lines are considered to have a 
Voigt profile (except for the Na and K doublet) and no trun¬ 
cation of the lines at large distances from the line cores has 
been applied. This choice of the far-wing treatment of the 
line shape is arbitrary. It is well known that the molecular 
lines should become sub-Lorentzian at large distances from 
the line core (see e.g., Freedman et al. 2008, and the refer¬ 
ences therein). However, the position of the cut-off, and the 
line wing shape itself, depend on the pressure and temper¬ 
ature and the perturber gases which broaden the molecular 
and atomic transitions. The choice to use Voigt profiles and 
not truncating the lines is thus only made because of the lack 
of knowledge regarding the actual line profiles. Grimm and 
Heng (2015) show that the differences when applying no line 
cut-off, when compared to an arbitrary cut-off, are at least of 
the order of 10 % when considering layer transmissions. In 
order to calculate the Voigt profiles we use the code provided 
by Humlfcek (1982). 

The calculations are performed on a pressure-temperature 
grid with 10 grid points in pressure going from 10“^ to 10^ 
bars (equidistantly spaced in log-space). Because the line 
wing strength due to pressure broadening is well behaved (the 
strength is simply linear in P), we found this grid spacing to be 
sufficient when interpolating to the actual pressures of inter¬ 
est. The temperature grid consists of 10 points going from 200 
to 3000 K, equidistantly spaced in log-space as well. Opac¬ 
ities with temperatures < 270 K are only calculated up to 1 
bar, temperatures up to 670 K only up to 10 bar, and tem¬ 
peratures up to 900 K only up to 100 bar. This choice was 
made because it was found, using the simple Guillot (2010) 
atmospheric model, that even cold planets such as Jupiter and 
Uranus should not be cooler than 270, 670 and 900 K at the 
pressures cited above. As we concentrate on hot jupiters in 
this paper we therefore did not extend the grid to cool tem¬ 
peratures at high pressures. We plan to extend the grid in the 
future, however. In total the above considerations yield 87 
pressure-temperature grid-points. 

Our fiducial wavelength range goes from 110 nm to 250yum. 
We calculate the opacities in this range on a grid with a spac¬ 
ing of A/AA =10^. This resolution is sufficient to resolve the 
line cores at all pressures and temperatures. From these calcu¬ 
lations we construct opacity distribution tables (k-tables) for 
later use (see Section 2.3.1). These tables are then interpo¬ 
lated to the pressure-temperature values of interest. 

Most of the line lists are obtained from the HI- 
TRANjHITEMP databases (Rothman et al. 2013, 2010), to¬ 
gether with additional data from the VALD, Kurucz and Ex- 
oMol line lists (Piskunov et al. 1995; Kurucz 1993; Harris 
et al. 2006; Barber et al. 2014). For methane we use HI- 


TRAN for temperatures below 300 K. For temperatures above 
300 K the ExoMol cross-sections are used (Yurchenko and 
Tennyson 2014), as this line list is much more complete at 
higher temperatures. The ExoMol cross-sections (in units of 
cm^ molecule"^) can be obtained, tabulated as a function of 
wavelength, from the ExoMol website. No pressure broad¬ 
ening has been applied when calculating these cross-sections, 
as pressure broadening information is not readily available. 
However, due to the sheer number of methane lines the cross- 
sections should be dominated by the Gaussian line cores. In 
general, for all molecular and atomic line opacities, pres¬ 
sure broadening information is often not available, especially 
when taking into account arbitrary mixtures of various molec¬ 
ular and atomic gaseous species. We therefore estimate the 
pressure broadening by using the air broadening coefficients 
7air of the HITRAN/HITEMP database when these are avail¬ 
able for a given molecule of interest. In cases where this infor¬ 
mation is missing as well we resort to the use of the pressure 
broadening approximation provided by Eq. (15) in Sharp and 
Burrows (2007). 

A special line shape treatment is needed when considering 
the Na (589.16 & 589.76 nm) and K (766.7 & 770.11 nm) 
doublet lines. Na and K are very important to correctly de¬ 
scribe the atmospheric absorption in the optical, as these two 
species are one of the main absorbers in this spectral range 
and their line wings act as a pseudo-continuum contribution 
to the total opacity (see, e.g.. Sharp and Burrows 2007; Freed¬ 
man et al. 2008). Different groups have tried to estimate the 
line shapes for Na and K taking into account collisions with 
H 2 and He (Burrows and Volobuyev 2003; Allard et al. 2003; 
Zhu et al. 2006), and the efforts are ongoing (Allard et al. 
2012). In particular Allard et al. (2003) showed that for brown 
dwarfs the use of correct Na and K wing profiles improves the 
agreement between synthetic spectra and observations. The 
line profiles we use for Na were obtained from Nicole Al¬ 
lard (private communication) using Rossi and Pascale (1985) 
pseudo potentials. For K we use the profiles available on the 
website of Nicole Allard^, which include C 2 v and Coov in¬ 
teraction symmetries. As H 2 should be the main perturber 
for alkali atoms in the atmosphere of giant planets, only H 2 - 
broadening is currently considered. The other lines of Na and 
K, which are much weaker than the doublet transitions, are 
modeled using van der Waals (vdW) broadening as described 
in Schweitzer et al. (1996): 

Tvdw = (1) 

where Ce is the van der Waals interaction constant, v is the 
mean relative velocity between H 2 and the alkali atom, Nu 2 
is the H 2 number density and a is a dimensionless number. 
Schweitzer et al. (1996) report that a = 11, but it is a factor of 
10 smaller in Sharp and Burrows (2007). We found that if we 
want to reproduce the vdW line widths given in Allard et al. 
(2007), then we need to use the smaller a value. The required 
ionization energies were taken from the NIST database^. 

Because the non-Lorentzian line profile calculations for the 
Na and K wings by Nicole Allard are only valid up to a certain 
H 2 density (10^^ or 10^^ cm“^ for Na or K, respectively) we 
revert to the use of Voigt profiles for higher densities. This 
occurs in the range of ~ 3-30 bars in hot jupiters, where the 
atmosphere becomes optically thick in the IR and the stellar 
light has been absorbed. 

^ http://mygepi.obspm.fr/~allard/alkalitables.html 

^ http: //WWW .nist.gov/pml/data/asd.cfm 
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Figure 1. Molecular and atomic line opacities of our database at a temper¬ 
ature of 1650 K and a pressure of 100 bar in our fiducial wavelength range 
going from 110 nm to 250 jum. Every color stands for a different species, 
with the names of the species indicated in the plot. 

In Figure 1 we show all molecular and atomic line opacities 
of our database at a temperature of 1650 K and a pressure of 
100 bar in our fiducial wavelength space going from 110 nm 
to 250 jum. The pressure of 100 bars is far higher than where 
the radiation in the planetary emission spectra usually stems 
from. However, as the pressure broadening smoothes out in¬ 
dividual lines, the large scale opacity features can more easily 
be seen at higher pressures. This figure has been generated 
from our opacity distribution database and shows the mean 
value Ky = Yji of every wavelength subgrid which are 
spaced at a resolution of A/AA = 1000. Here, g is the cumula¬ 
tive opacity distribution function, see Section 2.3.1 for more 
information. 

VO and TiO opacities have not been added yet. We ex¬ 
plained in Section 1 that the role of these two absorbers is 
quite controversial, as they might not be present in the at¬ 
mospheres due to a potential rain-out, cold-trap capture or 
photodissociation. Nonetheless, we plan to add VO and TiO 
opacities in the next version of the code. 

2.1.2. Continuum opacities 

As a continuum opacity source we currently consider col¬ 
lision induced absorption (CIA) arising from H 2 -H 2 and H 2 - 
He collisions. Tabulated data and programs from Borysow 
et al. (1988, 1989); Borysow and Frommhold (1989); Bo¬ 
rysow et al. (2001); Borysow (2002) were used to obtain the 
cross-sections^. 

2 . 2 . Stellar spectra 

For the host star spectral templates we use PHOENIX mod¬ 
els of main-sequence stars which have evolved to 1/3 of their 
main sequence lifetime.^ For the stellar evolution we use 
Yonsei-Yale tracks (Yi et al. 2001; Kim et al. 2002; Yi et al. 
2003; Demarque et al. 2004) as well as the evolutionary cal- 

^ Tables and code were obtained from http://www.astro.ku.dk/ 
-aborysow/programs/index.html 

^ The results for the stellar spectra depend only very mildly on this choice, 
the main effect being that the stars slowly increase their luminosity with time. 
Because the transiting hot jupiters that can be best studied orbit K-type stars, 
which typically have ages less than half of their main sequence lifetime, we 
chose a value of 1/3. 


culations of Baraffe et al. (1998). More details can be found 
in van Boekel et al. (2012). 

2.3. Code structure and modules 

The basic principle for solving for the atmospheric struc¬ 
ture is based on Dullemond et al. (2002), which we adapted 
to the case of 1-d plane parallel planetary atmospheres. The 
code starts with an initial guess for the temperature struc¬ 
ture, computes the corresponding molecular and atomic abun¬ 
dances and the resulting opacities and then calculates the tem¬ 
perature assuming radiative-convective equilibrium. The code 
then starts again with the newly found PT- structure until the 
solution converges. Because the PT -structure, the abundances 
and opacities mutually depend on each other, we solve for at¬ 
mospheric structure in an iterative fashion. 

From a given atmospheric temperature structure we obtain 
the molecular abundances using the CEA equilibrium chem¬ 
istry code (Gordon and McBride 1994; McBride and Gordon 
1996). When the current opacities are calculated, we solve 
the full angle and frequency dependent radiative transfer prob¬ 
lem of the planetary radiation field. From this we obtain the 
intensity-mean, fiux-mean and Planck mean opacities as well 
as the variable Eddington factors. These opacities and Ed¬ 
dington factors are then used in the Variable Eddington Fac¬ 
tor (VEF) module to find the temperature structure using the 
moments of the radiation field (see, e.g., Hubeny and Mihalas 
2014). The temperature is found using a two-stream approx¬ 
imation for the planetary and stellar radiation field. Further¬ 
more we check if a given atmospheric layer is convective by 
applying the Schwarzschild criterion. We switch to an adia¬ 
batic temperature gradient in the layers that are found to be 
convective. 

The iteration is stopped once the maximum change in tem¬ 
perature between the current iteration and the temperature 
found 60 iterations ago is smaller than 0.01 K and if the plane¬ 
tary emerging flux obtained from the full angle and frequency 
dependent radiative transfer solutions is equal to the imposed 
total flux with a relative maximum deviation of 0.001. In 
rare situations the iteration will slowly oscillate within a a 
given temperature range and not find a solution and therefore 
not converge to a solution with a maximum flux deviation of 
0.001. In this case we flag the files of the atmospheric struc¬ 
tures with _nconvergence_YYY where YYY is the relative de¬ 
viation to the imposed flux in percent. 

In sections 2.3.1, 2.3.2 and 2.3.3 we give further detailed 
information on the code modules (full radiative transfer mod¬ 
ule, VEF temperature iteration module and the equilibrium 
chemistry module, respectively). 

2.3.1. Radiative transfer module 

In protoplanetary disks, for which the VEF based approach 
by Dullemond et al. (2002) has been developed, the main 
contribution to the total opacity is coming from dust and ice 
grains (see, e.g., Semenov et al. 2003). An important feature 
of dust grain opacities is that they vary only slowly with wave¬ 
length, making it possible to use a small number of wave¬ 
length grid points. 

In the case of planetary atmospheres molecules contribute 
strongly to the total opacity. The molecules which are impor¬ 
tant for the spectrum can have hundreds of millions to tens 
of billions of very sharply peaked lines (see, e.g., Rothman 
et al. 2010; Yurchenko and Tennyson 2014, for H 2 O and CH 4 , 
respectively). Therefore, especially at low pressures, radia- 
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live transfer calculations need to be carried out at high spec¬ 
tral resolution. We thus have adopted the opacity distribution 
method and the correlated-k (c-k) assumption (Goody et al. 
1989; Fu and Liou 1992; Lacis and Oinas 1991) to carry out 
our radiative transfer calculations. Opacity distribution tables 
(k-tables) should yield a good description of the detailed high 
resolution opacities while keeping the numerical costs of the 
radiative transfer calculations minimal. 

We combine the k-tables of all molecular species contribut¬ 
ing to the total opacity using a fast combination method which 
has a computational cost linear in the number of species A^sp. 
see Section B. 

The calculations to obtain the mean opacities and Edding¬ 
ton factors are carried out on our fiducial grid (going from 
110 nm to 250 yum) with a grid spacing of T/AT =10, which 
results in 78 spectral bins. In order to test the accuracy of 
these results we have carried out calculations at a grid spac¬ 
ing of T/A/l = 50 as well, but the differences in the results are 
negligible. 

To calculate the emission spectrum of an atmosphere after 
the PT -structure has converged we carry out a c-k radiative 
transfer calculation at a grid spacing of A/AT = 1000, which 
results in 7729 spectral bins. 

In every spectral bin of the T/AT = 10, 50 and 1000 cases 
we employ a g-grid. g replaces the spectral coordinate T or 
y in c-k, it is equal to the value of the cumulative opacity 
distribution function, i.e. 


dg = f{K)dK , 


( 2 ) 


where f{K)&K is the fraction of the opacity values between k 
and a: -I- dA: within a given frequency interval. We carry out 
the radiative transfer calculation using g, instead of the fre¬ 
quency or wavelength. The frequency averaged mean Qy of 
any radiative quantity Qy within the spectral bin can then be 
calculated as 




1 

A^ 

f 


yv 

Jv 


Qydv' 


Qgdg, 


(3) 


where Qg is the quantity corresponding to Qy in g-space 
within the spectral bin of interest. 

For the T/AT = 1000 case we approximate g on a grid of 20 
Gaussian quadrature points, consisting of a 10 point Gaussian 
quadrature grid going from g = 0 to g = 0.9 and a 10 point 
Gaussian quadrature grid going from g = 0.9 to g = 1. 

For the T/AT =10 and 50 cases we take a finer grid in g. 
The g-grid has 36 points, consisting of a 6 point Gaussian 
quadrature grid ranging from g = 0 to g = 0.95, an 8 point 
Gaussian grid ranging from g = 0.95 to g = 0.99, a 20 point 
Gaussian grid ranging from g = 0.99 to g = 0.99999 and a two 
point trapezoidal quadrature grid ranging from g = 0.99999 
to 1. 

The different methods for obtaining the combined c-k opac¬ 
ity of all species at the resolutions of T/AT = 1000, A/AA = 
10 and A/AA = 50 are explained in Section B. 

The radiative transfer calculations are made using a 2nd or¬ 
der Feautrier method, considering 20 ju (i.e. cos 6) angles on 
a 20-point Gaussian quadrature grid. 


2.3.2. Variable Eddington factor method 



- A/AA =1000 — A/AA =50 . A/AA =10 

Figure 2. Upper panel: Emission flux of a hot jupiter calculated with the 
code introduced in this paper. The gray solid line shows the full line-by- 
line radiative transfer calculation at a resolution of T/A/l = 10^. Overplotted 
one can see the correlated-k calculations at T/A/l =10^ (black dashed line), 
T/A/l = 50 (red long dashed line) and at AIAA = 10 (blue short dashed line). 
Lower panel: Relative error of the AjAA = 10 ^, 50 , 10 calculations when 
comparing to the rebinned A!AA =10^ calculation. 

A description of the variable Eddington factor method, and 
how we use it to find the temperature in the atmosphere, can 
be found in Section C. 

In our code the radiation from the star can be received in 3 
different ways: (i) The angle between the atmospheric vertical 
and the stellar irradiation is (ii) The stellar flux 

is absorbed by the planet with a cross-section of but dis¬ 
tributed over the day side hemisphere (day side average): The 
incident vertical irradiation is reduced by a factor of 1/2; (iii) 
The stellar flux is absorbed by the planet with a cross-section 
of 7r/?pj but distributed over the full area (global aver¬ 

age): The incident vertical irradiation is reduced by a factor 
of 1 /4. In the day side or global average cases the stellar irra¬ 
diation field is treated to be shining at the atmosphere isotrop¬ 
ically. For our atmospheric grid we chose option (ii). We 
therefore assume that there is no efficient redistribution of the 
insolation energy to the night side. We will revisit the validity 
of this assumption in Section 3. 

2.3.3. Equilibrium Chemistry 

We use the NASA Chemical Equilibrium with Applications 
{CEA) code by Gordon and McBride (1994); McBride and 
Gordon (1996). The code minimizes the total Gibbs free en¬ 
ergy of all possible species while conserving the number of 
atoms of every atomic species. Given a pressure and tem¬ 
perature together with the atomic composition of the gas the 
output of the code is the mass and number fraction of all pos¬ 
sible outcome species (atoms, ions, molecules), the resulting 
density, as well as the adiabatic temperature gradient Vad of 
the gas mixture. 


2.4. Testing the code 

To characterize the quality of the results produced by the 
PETIT code a series of tests were carried out. 

2.4.1. Correlated-k radiative transfer 

First it was tested whether the correlated-k opacity combi¬ 
nation methods introduced in Section B yield results of suffi- 
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Figure 3. Bolometric flux of the converged atmospheric structure calculated 
from the d/A/l = 10 correlated-k radiation held integrated over angle and 
frequency space. The bolometric flux is shown as a black solid line. The 
two red solid vertical lines denote the imposed total and internal fluxes of the 
planet. The red shaded area denotes the radiative region of the atmosphere, 
whereas the blue shaded region shows the convective region. 


cient accuracy. To this end we calculated the emission spec¬ 
trum of a hot jupiter at our three different resolutions T/A/l = 
10^, 50, 10, using correlated-k and compared to the results 
of a line-by-line calculation at a resolution T/A/l = 10^. As 
an example PT -structure we took a self-consistent result from 
our code for a 1 M^, 1 planet^ around a sun-like star with 
an effective temperature = 5730 K with radius = R©. 
The planet was assumed to be in a circular orbit at a distance 
of J = 0.04 AU, have an internal temperature Tint = 200 K 
and a C/0 ratio of 1.17. We calculated the PT -structure for a 
day-side averaged hemisphere. 

The resulting emission spectra of the planet can be seen in 
the upper panel of Figure 2. In the lower panel we calcu¬ 
late the relative errors of the correlated-k calculations when 
compared to the frequency averaged line-by-line calculation. 
If the c-k assumption was perfectly valid the error would be 
zero, as the flux values of a c-k calculation at resolution (e.g.) 
10 should be identical to the flux of a higher resolution line- 
by-line calculation, after having been frequency averaged to 
the same resolution. One sees that in regions of appreciable 
flux the relative deviation between the rebinned 4/AT =10^ 
line-by-line calculation and the correlated-k calculations is al¬ 
ways smaller than 5 % and usually much less. Thus our results 
are within the accuracy limits commonly found for correlated- 
k (see, e.g., Fu and Liou 1992; Lacis and Oinas 1991). 


HD 189733b 



Figure 4. Secondary eclipse measurement of HD 189733b using data from 
Knutson et al. (2012); Charbonneau et al. (2008); Agol et al. (2010); Swain 
et al. (2010); Grillmair et al. (2008). Spitzer photometric and spectroscopic 
points are shown using black squares and crosses, respectively. The HST 
spectra is shown as black dots. Our model for HD 189733b at C/O = 0.4 and 
[Fe/H] = 1.4 is shown as a red line. For the photometric points we overplot 
boxcar-averaged red points obtained from our spectrum. 

Furthermore, deep within the atmosphere, but at lower pres¬ 
sure than the radiative-convective boundary Rconv, the radi¬ 
ation field only needs to carry the internal flux of the planet. 
The reason for this is that all the stellar flux has been absorbed. 
One thus finds that 

^deep(^ ^conv ) = ■ (5) 

Even further down the PT -structure will eventually become 
convective such that the radiative flux becomes negligible 
when compared to the convective flux. In Figure 3 one can see 
the result obtained from integrating the angle and frequency 
dependent radiation held of the T/AT =10 correlated-k struc¬ 
ture calculation. The radiation held was integrated to yield the 
bolometric flux in the atmosphere as a function of pressure. It 
can be seen that the surface flux converges to Fimposed- Fur¬ 
thermore, at approximately 3 bar, the stellar flux has been ab¬ 
sorbed and the radiative flux is equal to crT^^^. At even higher 
pressures (R ~ 70 bar) the atmosphere becomes convective 
and the flux transported by radiation starts to dwindle. The 
radiation held thus behaves as expected and the converged so¬ 
lution indeed fulfils the input parameters of the problem. The 
relative difference between the converged solution of the total 
emergent flux and the imposed flux was 0.08 %. 

2.4.3. Comparison to data: HD 189733b 


2.4.2. Energy balance 

As a next step we tested whether the converged solution is 
consistent with the input parameters. This was done by check¬ 
ing whether the final PT -structure, together with the molec¬ 
ular abundances and their corresponding opacities gives rise 
to the correct total emergent flux. For a day-side averaged 
PT -spectrum the total emergent flux should be 


imposed 


= (T 


p + 

^ int ^ 




(4) 


' Where and are Jupiter’s mass and radius, respectively. 


In order to get to get a qualitative impression of the compa¬ 
rability of our calculations with actual data we chose to look at 
HD 189733b, as it has quite a lot of available measurements. 
We used the following data: Spitzer IRAC photometry: 3.6 
and 4.5 pm (Knutson et al. 2012), 5.8 pm (Charbonneau et al. 
2008) and 8 pm (Agol et al. 2010), Spitzer IRS broadband 
at 16 pm , Spitzer MIPS at 24 pm (both Charbonneau et al. 
2008), HST NICMOS spectroscopy (Swain et al. 2010) and 
Spitzer IRS spectroscopy (Grillmair et al. 2008). For the stel¬ 
lar, planetary and orbital parameters we adopted T^ = 5040 
K, R, = 0.756 R©, Rpi = 1.138 R^ (Torres et al. 2008), Mpi 
= 1.137 (Butler et al. 2006; Agol et al. 2010) and d = 
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0.031 AU (Butler et al. 2006). The comparison between the 
data and our model for HD 189733b assuming C/0 = 0.4 and 
[Fe/H] = 1.4 can be seen in Figure 4. We chose these values as 
they are within the Bayesan regions with the highest credibil¬ 
ity identified by Benneke (2015) for this planet. One can see 
that between 5 and 8 yum our model fits the IRS spectroscopy 
quite well, including the water feature at 6.6 yum, while it is 
somewhat too high at larger wavelengths. While the IRAC 
points for wavelengths below 8 yum are not fitted very well, 
the IRAC photometry point at 8 yum, the IRS broadband and 
MIPS photometry points are all well fitted by our model. Fur¬ 
ther, the water absorption features between 1.5 and 2.5 yum in 
our model seem to correlate with the location of maxima and 
minima in the HST data. The depth of the absorption features 
is much bigger in the HST data, however, although some of 
the values in the HST spectra are negative, which is unphysi¬ 
cal and related to the observational process. We conclude that 
the comparison between observations and our model seems to 
already work quite well in certain parts of the spectra. Dedi¬ 
cated fitting studies might improve the results further. 

3. ATMOSPHERIC PROCESSES IN HOT JUPITERS 

As outlined in Section 1 our goal is to set up an atmospheric 
grid for hot jupiters. The range of effective temperatures we 
study for these objects is extending from 1000 K to 2500 K. 
We discuss important physical effects which govern the atmo¬ 
spheres of this class of planets below and assess how well the 
PETIT code is able to describe them. 

• Chemistry 

We are using a chemical equilibrium model for obtain¬ 
ing molecular and atomic abundances in our code and 
we assess the viability of this assumption below. As 
outlined before, the knowledge of these abundances is 
crucial to construct the atmospheric opacities. 

There are different regions in hot Jupiter atmospheres, 
in which different chemical assumptions are fulfilled. 

In the deep regions of the atmosphere temperatures and 
densities are high. Therefore the chemical reaction 
timescales are short. Here the chemistry is in equilib¬ 
rium, i.e. the system is in a state of minimal Gibbs free 
energy. By definition an equilibrium chemistry code 
will then suffice to obtain the molecular abundances. 

Further, there are two more important effects, which are 
often summarized in the term “non-equilibrium chem¬ 
istry”: 

In higher portions of the atmosphere the density is 
lower and the gas is often at lower temperatures: un¬ 
der these conditions vertical eddy diffusion can quench 
the abundances if the timescale for attaining chemical 
equilibrium is longer than the vertical mixing timescale. 

In even higher regions the density is very low. Here 
photodissociation, i.e. photochemistry can become the 
governing process if the insolation of the atmosphere 
is strong enough. In these regions the photodissocia¬ 
tion timescale will be shorter than the relevant chemical 
timescales. 

It is obvious that if the effective temperature (which 
translates into a distance to the star) is high, photodis¬ 
sociation acts on ever shorter timescales. However, this 
is compensated by the fact that a hotter atmosphere of 
a planet closer to its star will have shorter chemical 


timescales, such that planets at smaller semi-major axes 
are actually less affected by photochemistry than plan¬ 
ets further outside. 

We study how strongly quenching and photochemistry 
are expected to affect hot jupiters below. Our emphasis 
is on the regions which will shape the spectral appear¬ 
ance of the planet, i.e. the spectrally active regions. 

For emission spectra the spectrally active region of a 
planetary atmosphere usually lies in the pressure range 
from 10“^ to 10 bars (see, e.g., supplementary material 
in Madhusudhan et al. 2011). We obtain a reasonable 
assessment of the importance of non-equilibrium chem¬ 
istry by considering the work by Miguel and Kalteneg- 
ger (2014), who analyzed the chemical properties of 
planetary atmospheres around FGKM-stars. They used 
stellar model spectra compiled by Rugheimer et al. 
(2013) for the FGK stars and a spectral model for an 
inactive M dwarf by Allard et al. (2001). Furthermore 
they consider vertical mixing. 

By comparing figures 6 and 7 in Miguel and Kalteneg- 
ger (2014), we identify the effective temperature region 
where the spectrally active region is not affected by 
non-equilibrium effects to be Teff g [1500 K, 2600 K]. 
Given that only little stellar light is absorbed in the re¬ 
gions above 10“^ bar we do not expect the PT -structure 
for P > 10“^ bar to be compromised within this Teff- 
range.^ However, we want to remind the reader that 
our above choice of the Tgff-range is subjected to the 
assumptions made in Miguel and Kaltenegger (2014), 
in particular concerning the stellar model spectra, the 
eddy diffusion parameter (they took 10^ cm^ s"^) and 
the analytical model used for the atmospheric PT- 
structure. It is reassuring, however, that also Venot 
et al. (2015) did not find any significant differences 
in the emission spectra of hot C-rich planetary spec¬ 
tra when comparing different chemical schemes for the 
treatment of photochemistry. In their paper they com¬ 
pare a more sophisticated chemical network to the re¬ 
sults of a less complete carbon-chemistry network. Al¬ 
though the abundances of methane obtained by using 
these two different networks can vary by roughly an 
order of magnitude it leaves the emission spectra in 
their calculations unchanged (although methane is one 
of the strongest absorbers in these atmospheres). The 
reason for this is that the region where photochemistry 
becomes important lies above the spectrally active re¬ 
gion of the atmosphere. We thus conclude the atmo¬ 
spheres of the lowest temperatures in our grid could be 
affected by non-equilibrium chemistry. We thus flag the 
file names of all results with T^q < 1500 K with the 
“_neqc” flag to make the user aware of this. 

• Clouds 

Clouds appear to be widespread in all planetary at¬ 
mospheres. The most commonly stated evidence for 
clouds or hazes in hot Jupiter atmospheres is the fact 
that the transmission spectra of many of these ob¬ 
jects show no or only weak features at optical wave¬ 
lengths. This is strildng as in general one would ex¬ 
pect strong features from Na and K absorption in the 

^ We found that even in cases where the atmospheres are enriched in metals 
by up to 25 % (in mass) less than 10 % of the incident stellar radiation has 
been absorbed above the 10“^ bar altitude. 
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case of cloud free atmospheres. HD 189733b repre¬ 
sents a very prominent example, featuring a nearly flat 
transmission spectrum at optical wavelengths, except 
for the alkali line cores (e.g. Sing et al. 2011). Further 
(potential) examples for clouds or hazes weakening ab¬ 
sorption features in hot Jupiter transmission spectra are 
HD 209458b (Charbonneau et al. 2002), XO-2b (Sing 
et al. 2012), WASP-29b (Gibson et al. 2013a), HAT-P- 
32b (Gibson et al. 2013b) and WASP- 6 b (Jordan et al. 
2013). 

Clouds in hot jupiters may consist of silicates such as 
MgSiOs or Mg 2 Si 04 , liquid iron droplets, corundum 
(AI 2 O 3 ) and others. A further possibility is the photo¬ 
chemical creation of hydrocarbon hazes, arising from 
the photodissociation of CH 4 in the upper layers of the 
atmosphere. For a more detailed discussion of possible 
cloud and haze forming species see, e.g., Marley et al. 
(2013). 

Assessing the influence of clouds on the PT -structure 
and emission spectrum of hot jupiters is not an easy 
task. In the case of HD 189733b, which shows a fea¬ 
tureless optical transmission spectrum (except for the 
alkali line cores, see Sing et al. 2011), Barstow et al. 
(2014) And that the Pr-structure they can retrieve us¬ 
ing the planet’s emission spectrum is more or less in¬ 
sensitive to whether or not a cloud model is included 
(they use various MgSiOs models). At the same time 
many of their cloud models are able to reproduce HD 
189733b’s transmission spectrum. This indicates that 
for hot jupiters, at least for HD 189733b, the treatment 
of clouds is important for the appearance of the planet’s 
transmission spectrum, but not so much for the actual 
absorption of the bulk of the stellar light in the deeper 
layers of the day side atmosphere. In this case the in¬ 
fluence of clouds on the PP-structure and the emis¬ 
sion spectrum would be minor. This is in agreement 
with the earlier work by Fortney et al. (2008), who 
also And that clouds have a minor effect on their self- 
consistently calculated PT-proflles and emission spectra 
of hot jupiters and therefore neglect clouds. The obvi¬ 
ous importance of clouds in the case of transmission 
spectroscopy is due to the slant optical depths of possi¬ 
ble cloud species being -35-90 bigger than the vertical 
optical depth (Fortney 2005). 

We do not currently consider the formation of clouds 
and the associated effect on the planet’s opacity. How¬ 
ever, from the previous discussion we conclude that it 
might be permissible to neglect clouds in our calcula¬ 
tions. Nonetheless we want to note, following Fortney 
(2005), that in cases of high metallicity planets the ef¬ 
fects of clouds may become important, especially if ap¬ 
preciable amounts of silicate, iron or corundum con¬ 
densates can form. This has to be stressed in light of 
the fact that hot jupiters seem to be most prevalent in 
stellar systems of high metallicity (Fischer and Valenti 
2005). 

• Winds 

Based on GCM simulations and theoretical considera¬ 
tions, winds are expected to be present on hot jupiters, 
driven by the temperature contrasts between the day 
and nightside, and the polar and equatorial regions 
(see, e.g., Heng and Showman 2014). The question 


of whether these winds will have an effect on the ther¬ 
mal structure of the planetary atmosphere depends on 
whether the advection timescale of the winds Tadv is 
shorter than the radiative cooling timescale Trad and/or 
chemical timescale Tchem of the atmosphere. If Tadv is 
indeed shorter than one of those two timescales, then 
energy or molecules will be transported, and the as¬ 
sumptions of local radiative or chemical equilibrium 
breaks down. To properly carry out this timescale com¬ 
parison one would have to couple GCM simulations 
with radiative transport and chemical non-equilibrium 
calculations, which is beyond the scope of this work. 
The fact that one sees a day-night temperature varia¬ 
tion when looking at the thermal phase curve of, e.g., 
HD 189733b (Knutson et al. 2012), shows that winds 
are not able to perfectly redistribute the energy from 
the incident stellar radiation across the whole planetary 
surface. However, the results in Knutson et al. (2012) 
also show that the hottest and coldest points in the at¬ 
mosphere are offset from the substellar and antistellar 
point, respectively. This indicates that winds play a role 
in distributing energy across the planet. In general, it 
is found that the higher the effective temperature of a 
hot Jupiter, the less efficient the transport of energy by 
wind becomes (Perez-Becker and Showman 2013). For 
“cool” planets with effective temperatures of ~ 1000 K 
redistribution of energy may be quite efficient unless 
the planet has a mass of a few Jupiter masses or more 
(Kammer et al. 2015). 

In order to at least partially accommodate the effect of 
heat redistribution by winds, our code has 3 possible 
ways to treat the distribution of the incident stellar fight 
across the atmosphere: (i) no wind transport of energy, 
(ii) day-side averaging or (iii) global averaging, the lat¬ 
ter approximating the case where winds highly effi¬ 
ciently distribute the energy received by the star across 
the planetary surface (see Section 2.3.2). Our treatment 
of the stellar energy input in the cases (ii) and (iii) are 
only approximative ways to inject the stellar energy into 
the planetary atmosphere. A fourth way would be to 
use a redistribution parameter for the incident stellar ir¬ 
radiation which adds a fraction of the absorbed stellar 
energy to the night side internal temperature and de¬ 
creases the amount of fight to be absorbed on the day- 
side (Burrows et al. 2006). Other possibilities include 
the mimicking of planetary winds by assuming that the 
atmosphere carries out a rigid body rotation, as it was 
done in Iro et al. (2005). 

4. SETUP AND CALCULATION OL THE GRID 
4.1. Grid setup 

We set up a grid of 10,640 models which is defined by the 
following parameters: 

1. TefF = 1000,1250,1500,1750, 2000,2250, 2500 K 

We chose to go to temperatures somewhat lower than 
where we are unaffected by non-equilibrium chemistry 
effects (1500 K, see Section 3). The files of models with 
TefF < 1500 K will be flagged with “_neqc” to make the 
user aware of potential differences when including non- 
equilibrium chemistry. Furthermore, high metallicity 
models with low log(g) and high TgfF will have temper¬ 
atures larger than 3000 K in the higher pressure parts 
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of the atmosphere. If this happens before the atmo¬ 
sphere becomes convective we flag these models with 
“_t3000k”, as our opacity grid only extends to 3000 
K (see Section 2.1.1). At atmospheric layers where 
T > 3000 K we use the opacities at 3000 K. 

2. [Fe/H] = -0.5, 0.0, 0.5,1.0, 2.0 

The metallicity is chosen to range from slightly subso¬ 
lar to strongly enriched and we use scaled solar compo¬ 
sitions according to Asplund et al. (2009). It is not gen¬ 
erally expected that enriched exoplanets have a scaled 
solar composition. Nonetheless, we use this approxi¬ 
mation as a proxy for various degrees of enrichment. 
A further degree of freedom regarding the composition 
is introduced to our grid by varying the C/0 ratio. In 
this work we focus on metallicities higher than the so¬ 
lar value. The reason for this is that giant exoplanets 
are expected to be enriched in metals, with objects of 
several hundred Earth masses having metallicities of up 
to several tens of the solar metallicity (Fortney et al. 
2013). 

3. C/O = 0.35,0.55,0.7,0.71,0.72,0.73,0.74,0.75,0.85, 
0.9,0.91, 0.92, 0.93, 0.94,0.95,1.0,1.05,1.12,1.4 

We investigate C/O values which are subsolar or su¬ 
persolar but < 1 (C/O© ~ 0.55), as well as values 
around and above 1. We use a flner sampling around 
C/O ~ 0.73 and C/O ~ 0.92, because we want to re¬ 
solve the transition from oxygen to carbon-dominated 
spectra and atmospheres at low and high temperatures. 
Commonly, the transition is expected to happen quite 
sharply at C/O values around 1 (see, e.g., Kopparapu 
et al. 2012; Madhusudhan 2012). We And C/O = 0.92 
for the high temperature atmospheres. Furthermore, the 
infrared opacity of the atmospheres is minimal when 
C/O is close to 1, because most of the C and O atoms 
are locked up in CO and neither H 2 O nor CH 4 of HCN 
are very abundant. This gives rise to inversions for the 
hottest atmospheres (Tgff > 1500) K, where the alkali 
atoms absorb the stellar irradiation quite effectively but 
the cooling is inefficient due to the IR opacity mini¬ 
mum (see Section 5). The C/O ratio at a given metallic¬ 
ity was obtained from varying the O abundance. This 
means that for supersolar C/O ratios the O abundance 
was decreased, corresponding to the accretion of water 
depleted gas or planetesimals during the planet’s for¬ 
mation. 

4. Spectral type of host star: F5, G5, K5, M5 

In order to assess the dependence of the atmospheric 
structure on the spectral shape of the stellar radiation 
field we calculated our grid using 4 different spec¬ 
tral types for the host star. For the earlier spectral 
types the energy received by the planet is absorbed pre¬ 
dominantly by the alkalis in the optical wavelengths, 
whereas for the later spectral types the wavelength 
range of the absorption shifts more and more to the 
IR, leading to increasingly isothermal planetary atmo¬ 
spheres. 

5. log(g) = 2.3,3.0, 4.0,5.0 

Our log(g) grid was chosen such that it encompasses 
hot jupiters of every conceivable mass-radius combina¬ 
tion, including bloated hot jupiters as well as compact 
(Rpi ~ R^) planets of varying masses (all planets listed 


on http: //exoplanets. org with a mass and radius 
measurement fall within our adopted log(g) range). 

4.2. Chemical model 

The following atomic species were considered in the equi¬ 
librium chemistry network: H, He, C, N, O, Na, Mg, Al, Si, P, 
S, Cl, K, Ca, Ti, V, Fe and Ni. Based on Seager et al. (2000) 
we consider the following reaction products: e“, H, He, C, N, 
O, Na, Mg, Al, Si, P, S, K, Ca, Ti, Fe, Ni, H 2 , CO, OH, SH, 
N 2 , O 2 , SiO, TiO, SiS, H 2 O, C 2 , CH, CN, CS, SiC, NH, SiH, 
NO, SN, SiN, SO, S 2 , C 2 H, HCN, C 2 H 2 , CH 4 , AlH, AlOH, 
AI 2 O, CaOH, MgH, MgOH, VO, VO 2 , PH 3 , CO 2 , Ti 02 , Si 2 C, 
Si 02 , FeO, NH 2 , NH 3 , CH 2 , CH 3 , H 2 S, KOH, NaOH, NaCl, 
KCl, H^, H“, Na^, K^, Fe (condensed), AI 2 O 3 (condensed), 
MgSi 03 (condensed), SiC (condensed). 

The choice of condensed species is motivated by Seager 
et al. (2000); Sudarsky et al. (2003). Additionally, we also 
added SiC as a condensable species, to account for conden¬ 
sation of C in atmospheres with a high C/O ratio, as has also 
been suggested by Seager et al. (2005). 

A reaction pathway that is of prime interest is the one con¬ 
necting H 2 O, CH 4 and CO. In Section 1 we already intro¬ 
duced the C/O ratio as a useful quantity for characterizing 
planetary atmospheres, as it allows to interpret the relative 
abundances of CH 4 and H 2 O for temperatures T < 1750 K. 
CH 4 and H 2 O are important molecules because they are abun¬ 
dant, have a high infrared opacity and therefore shape the 
overall appearance of the atmosphere’s emission spectrum. 

The net chemical equation of interest for this case is 

T > 1000 K 

CH4 + H2O . CO -r 3H2 , (6) 

T < 1000 K 

leading to a quite sharp transition of CH 4 vs. H 2 O rich atmo¬ 
spheres at C/O ~ 1 (see, e.g., Kopparapu et al. 2012; Mad¬ 
husudhan 2012 ): 

In chemical equilibrium CO is the most common C and O 
bearing molecule in planetary atmospheres, where the tem¬ 
peratures are high enough (T > 1000 K). In an oxygen-rich 
atmosphere (C/O < 1) the remaining oxygen is then partly 
found in the form of H 2 O and almost no CH 4 is present, as 
most C is locked up in CO. In a carbon-rich atmosphere (C/O 
>1) the excess C is put partly into CH 4 , with no O left to form 
water. For T < 1000 K the low temperature direction in Eq. 
( 6 ) is dominant, leading to appreciable amounts of both CH 4 
and H 2 O and negligible amounts of CO. 

The main effect of including condensation is the removal 
of oxygen from the gas phase through the condensation of 
MgSi 03 for temperatures smaller than ~ 1500 K, leading to a 
spectrally noticeable decrease of H 2 O and CO at C/O values 
as low as 0.7. This effect is observed for the Tgff = 1000, 
1250 and 1500 K cases. Effectively it shifts the spectrally 
visible transition from H 2 O dominated atmospheres to CH 4 
dominated atmospheres away from C/O ~ 0.9 to somewhat 
smaller values of C/O ~ 0.7, as the formation of MgSi 03 acts 
as a sink for the O atoms available to form H 2 O. 

The depletion of 0-bearing gas phase species due to con¬ 
densable 0 -bearing species has been found in much more 
complete cloud models as well (Helling et al. 2014). We de¬ 
scribe some of the incompletenesses of our cloud model be¬ 
low: One of the effects our condensation model does not treat 
is the the problem of homogeneous or heterogeneous nucle- 
ation, which could potentially shift the formation of conden¬ 
sates in the atmospheres toward layers of higher supersatura- 
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tion if initial condensation seeds are not present in the atmo¬ 
spheres (see, e.g., Marley et al. 2013). 

Further we want to stress that the condensed species in each 
layer remain in chemical contact with the gas phase in our 
model and do not rain out to deeper layers of the atmosphere. 

The consequences of a potential rainout for a planetary at¬ 
mosphere can be manyfold. First of all the rainout removes 
metals from the atmosphere, relocating them to deeper lay¬ 
ers. Hence the corresponding grain or droplet opacity will be 
missing from higher atmospheric layers. Because we do not 
include cloud opacities in our calculations we make the im¬ 
plicit assumption of a rainout of the condensed particles, al¬ 
though we do not model it, the net effect being the removal of 
metals from the higher layers. It has to be kept in mind, how¬ 
ever, that the chemical equilibrium solution of the gas abun¬ 
dances in chemical contact with the condensed species is not 
necessarily the same as it would be when assuming a rain¬ 
out. Our implicit assumption of a rainout is also applicable 
when considering the gaseous Na and K alkali abundances. 
In our models MgSi 03 condenses at temperatures below ~ 
1600 K. In principle this silicate material can further react 
with the alkali atoms to form alkali feldspars (such as albite 
and orthoclase), removing the gaseous alkalis from the gas 
for T < 1600 K (see, e.g.. Fodders 2010). We do not con¬ 
sider these feldspars in our condensation model, such that the 
alkali atoms stay in the gas, as they would in a silicate rain¬ 
out scenario. It has been found that alkali atoms are present in 
cool brown dwarf atmospheres, indicating that silicate rainout 
may occur in these objects (Marley et al. 2002; Morley et al. 
2012). Another consequence of condensed material can be 
the formation of a cloud deck, close to and above the layers of 
the atmosphere hot enough the evaporate the in-falling cloud 
particles again. Such cloud decks can heat the atmosphere lo¬ 
cally and in the layers below, by making the atmosphere more 
opaque to the planet’s intrinsic flux, effectively acting like a 
blanket covering the lower layers of the atmosphere (see, e.g., 
Morley et al. 2014; Helling and Casewell 2014). If the cloud 
layer is optically thick close to the planet’s photosphere it will 
leave an imprint on the planet’s spectral appearance and and 
may reduce the contrast of absorption features. The height of 
the cloud deck depends critically on the planets effective tem¬ 
perature and also on its surface gravity since the condensa¬ 
tion temperature is pressure-dependent. The cooler an object 
is, the deeper in its interior the clouds will reside. There¬ 
fore the spectral imprint of clouds will vary with temperature, 
similar to the behavior in brown dwarf atmospheres. Silicate 
clouds with a high optical depth are thought to reside in the 
photospheres L4-L6 type brown dwarfs (Tgff ~ 1500-1700 K) 
where they affect the spectra. For cooler objects the cloud 
deck lies below the photosphere and the clouds are no longer 
seen (see, e.g. Fodders and Fegley 2006). In our atmospheres 
we checked the possible locations of the cloud decks (i.e. the 
layers below which the condensates evaporate). We found that 
the silicate evaporation layer of planets with Tgff = 1000 K and 
Teff =1250 K is always located at pressures far higher than that 
of the photosphere, such that we do not expect any spectral 
impact of a cloud layer. For effective temperatures between 
1500 K and 1750 K the evaporation layer lies close to and 
above the photosphere (in altitude), such that a cloud deck 
could potentially affect the spectrum. For increasing log(g) 
the photosphere shifts to layers of deeper pressure, but so does 
the evaporation layer, as condensation is pressure dependent. 
Note that this temperature range is close to the effective tem¬ 
perature where F4-F6 dwarfs are thought to be most strongly 


affected by silicate clouds. For higher temperatures the evap¬ 
oration layer is far above the photosphere such that we do not 
expect clouds to be of importance. 

For C/0 ratios > 1 and temperatures > 1750 K we find, in 
agreement with previous studies, that the spectrally most im¬ 
portant carbon bearing molecule is no longer CH 4 , but HCN 
(see, e.g., Kopparapu et al. 2012; Venot et al. 2012; Moses 
et al. 2013). In general the lower the pressure and the higher 
the temperature the more important HCN becomes. Therefore 
we see that the spectra at the highest effective temperatures are 
dominated by HCN absorption. 

4.3. Calculation of the grid 

The calculations were carried out using 150 atmospheric 
layers spaced equidistantly in log(R) between 10 “^"^ and 
9x10"^ bar. Note that our opacity grid is only calculated be¬ 
tween 10“^ and 10^ bar. For pressures outside this range we 
use the opacities at the pressures at the boundaries of our 
opacity grid. The grid calculations were extended to smaller 
pressures to not introduce any kinks at the 10 “^ boundary: the 
alkali line cores are already optically thick at these low pres¬ 
sures, and a cut-off of the atmospheric structure at 10 “^ bar 
would result in no alkali core flux coming from above at the 
highest point in the atmosphere, making the temperature there 
to cool. We provide the PT -structures only between 10“^ and 
10^ bar. However, at altitudes above the 10“^ bar level the 
contribution of the pressure-broadened line wings is to the to¬ 
tal opacity is negligible and the opacity is dominated by the 
line cores, whose shape is given by thermal broadening and 
is independent of pressure. As only little mass is above any 
given pressure lower than 10 “^ bar, the line wings are not able 
to significantly alter the radiation field. Therefore, adapting 
the 10 “^ opacity curves at all lower pressures should not af¬ 
fect the resulting PT structures; in all this range the line cores 
are of significant optical thickness, whereas the line wings are 
highly optically thin. Hence, the line cores govern both the 
absorption and the re-emission of energy, and thus the PT- 
structure. 

The calculations were extended to pressures larger than 10^ 
bar as we consider quite large surface gravities, which essen¬ 
tially rescale the temperature structures to higher pressures. 
We wanted to make sure that we do not cut-off the atmo¬ 
spheric structures at 10 ^ bar for high log(g) cases when the 
atmosphere is not yet optically thick at all wavelengths. We 
found, however, no differences in the PT structures nor the 
emission spectra when comparing cases extending down to 
either 10 ^ or 9x10"^ bar. 

For the temperature iteration the pressure-, temperature- 
and abundance-dependent combination of the individual 
species’ opacity tables is the computationally most demand¬ 
ing part of the atmospheric structure calculation. Thus, for nu¬ 
merical convenience, we precalculated the opacity tables for 
every atmospheric structure on 40 x 40 pressure and tempera¬ 
ture grid points (taking about 2 minutes) before the iterations 
were run. We then interpolated in this table during the itera¬ 
tions and verified that the results were consistent with those 
obtained when re-calculating the opacity tables for every in¬ 
dividual iteration. 

4.3.1. Convection and convergence 

As described in Section (C.3), we use the Schwarzschild 
criterion to assess whether a given layer in the atmosphere 
should be convective, and if so we switch to an adiabatic tern- 
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Figure 5. Atmospheric PT-structures for planets of varying host star spectral types, effective temperatures and C/O ratios with log(g) = 3 and [Fe/H] = 1. The 
line style varies with host star spectral type as follows: F5 (solid), G5 (dashed), K5 (dot-dashed), M5 (dotted). The line color indicates the following planetary 
effective temperatures: 1000 K (black), 1250 K (blue), 1750 K (purple), 2250 K (red). The four different panels correspond to 4 different C/O ratios: C/O = 0.55 
{upper left panel), C/O = 0.85 {upper right panel), C/O = 0.95 {lower left panel), C/O = 1.4 {lower right panel). 


perature gradient. We find that the lowest layers of the at¬ 
mospheres (at the highest pressure) become convective, with 
a radiative gradient much bigger than the adiabatic tempera¬ 
ture gradient. For hot atmospheres (Teff > 2000 K) with high 
metallicities [Fe/H] > 1 we find that regions with a steep tem¬ 
perature gradient high in the atmosphere (10“^ bar > P > 10“^ 
bar) can become convective. In these situations the solutions 
can become unstable, as the layers switch back and forth be¬ 
tween being either radiative or convective, introducing jumps 
and kinks in the PT-spectra. This suggests that these lay¬ 
ers are in the continous transition region between being fully 
radiative or convective, which cannot be resolved by the bi¬ 
nary Schwarzschild criterion. A better treatment would be to 
implement convection via the mixing length theory (MLT), 
as it allows for a continuous transition from a fully radiative 
to a fully convective solution. For now, we decided to re¬ 
run the PT -structures affected by this convergence problem 
and to forbid the occurrence of convection in the uppermost 
layers (10“^ bar > P > 10“^ bar) of the atmosphere. The cor¬ 
responding atmospheric structure files have been flagged with 
“_conv”. We plan to implement MLT in a future version of 
the code. 


5. RESULTS 

We first discuss some general characteristics of our results 
in Section 5.1. We will study the atmospheric properties sys¬ 
tematically as a function of effective temperature for all atmo¬ 
spheric parameters in sections 5.5 - 5.7. 


5.1. A first glance 

To give a first overview of our of results we show atmo¬ 
spheric PT -structures of log(g) = 3 and [Fe/H] = 1 planets 
for varying host star spectral types (F5, G5, K5, M5) and ef¬ 
fective temperature (1000 K, 1250 K, 1750 K, 2250 K) at four 
different C/O ratios (0.55, 0.85, 0.95, 1.4) in Figure 5. 

Some general, expected trends can quite easily be made out 
from looking at this plot: 

• The later the host star spectral type, the more isothermal 
the atmospheric structure becomes. This is expected 
because the wavelength range of the received stellar ir¬ 
radiation becomes more and more similar to the wave¬ 
length range of the internal planetary radiation field, 
such that the radiation field absorbed by the gas at the 
top of the atmosphere is similar to the radiation field 
absorbed by the gas at the bottom of the atmosphere. 
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Figure 6. Left panel: Pr-structure of a log(g) = 3, [Fe/H] = 1, Peff = 2250 
K, C/O = 0.95 atmosphere of a planet in orbit around a G5 star. Right panels: 
Local stellar flux (red solid line) at the three pressure levels at 3.47 x 10“^ 
(top panel), 9.07 x 10“^ (middle panel) and 1.27 x 10“^ bar (bottom panel) 
in the atmosphere. The local opacity log(/c:) for each layer is shown as a grey 
solid line (rescaled). The respective pressure levels are indicated by a red 
circle, square and diamond in the PT-structure. 

hence leading to similar temperatures. 

• The PT -structures with C/O = 0.55 are hotter than the 
RT-structures with C/O = 0.85. The main reason for 
this is that the atmosphere with the lower C/O ratio has, 
everything else being equal, more oxygen and thus a 
higher opacity due to a higher H 2 O abundance. This 
results in a stronger green house effect, as the excess 
H 2 O leads to a less efficient escape of radiation from 
the atmosphere. In order to radiate away the required 
amount of energy (set by Teff) the atmospheres need to 
be hotter. 

Another very striking result is that for C/O ratios close to 1 
temperature inversions form in the atmospheres for effective 
temperatures above 2000 K. In general, they can even occur at 
effective temperatures as low as 1500 K, see Section 5.6. This 
is interesting, as no extra optical opacity sources such as TiO 
and VO except for the ones given in Table 1 are being con¬ 
sidered. For host stars later than K5 there are no inversions in 
the planetary atmospheres. This phenomenon will be further 
studied in Section 5.1.1. 

5.1.1. Inversions at high CjO ratios 

As outlined above, C/O ratios of ~ 1 can lead to inversions 
in atmospheres with high enough effective temperature if the 
stellar host is of K spectral type or earlier. The reason for the 
inversions to occur for these spectral types is that an apprecia¬ 
ble amount of stellar flux is received from the star in the op¬ 
tical wavelength regime. This means that the alkali lines, and 
the pseudo-continuum contribution of the alkali line wings, 
will become very effective in absorbing the stellar irradiation. 

At the same time, close to C/O = 1, most of the oxygen and 
carbon is locked up in CO, leading to low H 2 O, CH 4 and HCN 
abundances and opacities. 

The combined effect of the effective absorption of the 
strong irradiation and a decreased ability of the atmospheric 
gas to cool, because of too little CH 4 , H 2 O and HCN leads to 
the inversion in the atmospheres. 

The absorption of the stellar light as a function of depth can 
be seen in Figure 6 , where we plot the PT -structure of a log(g) 
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Figure 7. Mass fractions of CH 4 (thin solid red line), H 2 O (dashed blue line), 
HCN (dotted green line) and CO (thick solid black line) as a function of the 
C/O ratio for a log(g) = 3, [Fe/H] = 1 atmosphere of a planet in orbit around a 
G5 star at a pressure level of 9.07 x 10“^ bar. The top panel shows the mass 
fractions for a planet with Teff = 1750 K while the bottom panel shows the 
mass fractions for a planet with Teff = 2250 K. 

= 3, [Fe/H] = 1, Teff = 2250 K, C/O = 0.95 atmosphere of a 
planet in orbit around a G5 star, as well as the local stellar 
flux at the pressure levels 3.47 x 10“^, 9.07 x 10“^ and 1.27 
X 10“^ bar in the atmosphere. Also a plot of the logarithm of 
the (rescaled) opacity log(A:) is shown in the Figure for each 
pressure level. The respective pressure levels are indicated by 
red points in the PT -structure. 

Figure 6 nicely shows how the alkali pseudo-continuum ab¬ 
sorbs the full stellar flux in its wavelength domain at the po¬ 
sition of the inversion: At the highest pressure shown in the 
spectral plots (3.47 x 10“^ bar) the stellar flux is still com¬ 
pletely unaffected by any absorption effects as the atmosphere 
is still optically thin at all wavelengths (except for right at the 
core of the alkali lines). At the hottest point in the tempera¬ 
ture inversion (at 9.07 x 10“^ bar) one can see that the alkali 
wings have already started to absorb non-negligible amounts 
of energy, and just after the inversion (at 1.27 x 10“^ bar) the 
stellar flux in the alkali wings has been completely absorbed. 
Interestingly, the inversions obtained in our calculations due 
to alkali heating seem to abide by the rule that the tropopause, 
i.e. the atmospheric layer at minimum temperature just after 
the inversion, should commonly be found at ~ 0.1 bar for a 
wide variety of possible atmospheres (Robinson and Catling 
2014). 

As can be seen in the stellar flux spectrum at the highest 
pressure the absorption of the stellar light outside of the alkali 
wings is rather sluggish, showing the importance of the alkali 
wings in the formation of the inversion. 

As mentioned above, in a small region of C/O around 1, the 
atmosphere’s ability to efficiently radiate away the absorbed 
stellar light decreases due to the involved chemistry. This can 
be understood by looking at Figure 7, which shows the CH 4 , 
H 2 O, HCN and CO mass fractions in a log(g) = 3, [Fe/H] = 1 
atmosphere of a planet in orbit around a G5 star as a function 
of C/O at a pressure level of 9.07 x 10“^ bar, i.e. close to 
the pressure where the inversion temperature, if an inversion 
occurs, is maximal. Two cases for planets with Teff = 1750 
K and Teff = 2250 K are shown and we carried out 100 self- 
consistent atmospheric calculations for both cases with C/O 
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Figure 8. Emission spectra as a function of host star spectral type for a Teff 
= 1750 K, log(^) = 3, [Fe/H] = 1 planet with C/O = 0.55 (upper panel) and 
C/O = 1.05 (lower panel). The spectra are shown for a F5 (blue lines), G5 
(purple lines), K5 (red lines) and M5 (orange lines) host star. The colored 
bars indicate the position of the absorption maxima of various species. The 
black line shows the blackbody flux at the atmosphere’s effective temperature. 

going from 0.35 to 1.4 in equidistant steps. 

One sees that for the Teff = 2250 K case, at C/O = 0.95, the 
H 2 O abundance has already decreased by 4 orders of magni¬ 
tude when compared to the lowest C/O values, while the CH 4 
abundance is still more than 2 magnitudes smaller than its 
highest abundance at the highest C/O values. Further, HCN 
has not yet risen to a high enough abundance to take over 
the cooling. The C/O = 0.95 point at Teff = 2250 K thus is 
very close to the aforementioned point of minimum IR opac¬ 
ity, leading to the inversions seen in our results for all host 
spectral types except M5. For higher C/O ratios the IR opac¬ 
ity and the atmosphere’s ability to cool increases, such that no 
inversions are observed anymore, mainly because HCN takes 
over the cooling. 

For the particular case of Teff = 1750 K in Figure 5 the situa¬ 
tion must be different, as there is no inversion present in the at¬ 
mosphere. The reason for this can be seen in the panel for Tgff 
= 1750 K in Figure 7: for this atmosphere the transition from 
water-rich to methane-rich atmospheres occurs much quicker 
as a function of C/O than it does for the Teff = 2250 case. The 
methane mass fraction jumps from 10“^ to 10“^ at C/O = 0.93 
and the HCN mass fraction jumps from 10“^ to 10“"^ and no 
extended region of low water, methane and HCN abundance 
is seen. Further, as this atmosphere is cooler, the overall CH 4 
content is higher than in the hotter case. This is expected to 
occur and has been studied before both in equilibrium and dis¬ 
equilibrium chemical networks (see, e.g., Moses et al. 2013), 
showing that CH 4 becomes less abundant as the temperature 
increases in carbon-rich atmospheres. In conclusion, this at¬ 
mosphere can cool more efficiently. 


5.1.2. Inversions and line list completeness for HCN and C 2 H 2 

We want to issue a word of caution regarding the cooling 
efficiency of atmospheres. At high temperatures for C/O > 
1 and Teif > 1750 K we find that HCN is more abundant 
than CH 4 . It is therefore very important to use HCN line lists 
which are as complete as possible. In fact we found that if we 
use HCN from the HITRAN database, which is made for low 
atmospheric temperatures, we got strong inversions occurring 
even for C/O > 1 if the effective temperatures were high. Only 
once we switched to the ExoMol line list for HCN we got the 


results presented in this paper, where inversions only occur 
for C/O ~ 1. The ExoMol line list is much more complete 
for HCN, containing many more lines. The line list is made 
specifically for high temperatures, optimized for temperatures 
up to 3000 K and compares well to a high temperature labora¬ 
tory measurement made at T = 1370 K (Barber et al. 2014).^^ 
This allows the atmospheres to cool more efficiently, making 
the inversions go away in many cases. 

Likewise, we want to stress that we use the HITRAN line 
list for the C 2 H 2 molecule, as an ExoMol version is not avail¬ 
able. C 2 H 2 is quite common in our results for C/O > 1 in the 
cases where HCN is common as well. This suggests that the 
atmospheres ability to cool might be further enhanced if high 
temperature line lists for C 2 H 2 were to be considered. 

5.2. Host star dependance of the atmospheres 
5.2.1. Spectra 

As described in Section 5.1 planets orbiting increasingly 
cooler host stars will approach an increasingly isothermal at¬ 
mospheric structure, because the spectral energy distribution 
of the insolation becomes more and more comparable to the 
SED of the planetary radiation field. 

We show the emission spectra of atmospheres with varying 
host star spectral type for a planet with Tq^ = 1750 K, log(g) 
= 3, [Ee/H] = 1 for two different C/O ratios (0.55, 1.05) in 
Eigure 8 . We indicate the positions of absorption features of 
H 2 O, CO 2 , K, Na, CO, CH 4 , PH 3 and HCN in the plots. Eor 
the atmospheres with C/O = 0.55 the emission spectra clearly 
become more blackbody-like as the host star gets cooler: the 
excess emission (with respect to the blackbody curve at 1750 
K) of the atmospheres for T < 1.3 pm decreases for cooler 
host stars. Eurthermore the molecular absorption bands in the 
emission spectra start to get shallower. As expected for a C/O 
ratio < 1 , the spectra are clearly water-dominated. 

Eor the atmospheres with C/O = 1.05 the situation is some¬ 
what different. Eirst, the atmospheres are clearly carbon- 
dominated, showing strong HCN features. Moreover, the lat¬ 
est type host star (M5) causes the least isothermal planetary 
spectrum, while all earlier type host stars result in a much 
more isothermal atmospheric structure and, therefore, spec¬ 
tra. This is the contrary of what we saw for the C/O = 0.55 
atmosphere, now host stars of an earlier type are making the 
planetary spectra more isothermal. This is merely the spec¬ 
tral consequence of early type host stars creating inversions 
or isothermal atmospheres for planets with C/O ~ 1, which 
we explained in Section 5.1.1. As the M5 star is not able to 
heat the atmosphere enough due to a lack of energy in the op¬ 
tical wavelengths the corresponding PT -structure and spectra 
are less isothermal. The FT-structures producing the spec¬ 
tra shown here for C/O = 1.05 do not have inversions, they 
are just more isothermal due to the heating. As we will see 
in Section 5.7, atmospheres at Teff = 1750 K can, in general, 
exhibit inversions. 

5.3. log(g) dependence of the atmospheres 
5.3.1. PT-structures 

The behavior of the PT -structures as a function of log(g) 
is studied in Eigure 9. If one considers gray opacities which 
are constant as a function of P and T and assumes hydrostatic 

More comparisons are not possible as there are not many high tempera¬ 
ture measurements for this molecule. 








14 


R Molliere et al. 



, , , I , , , , I V. , /l*", T ^ , I , T", . I .... I XU i .I . I .I .I . 

1000 2000 3000 4000 iq-^ iq^ iq’^ 10'^ 10'^ 10° 
Temperature (K) Stellar flux absorption 


Figure 9. Atmospheric PT -structures and stellar light absorption as a func¬ 
tion of log(g) for planets with Teff = 1250 K, [Fe/H] = 1 and C/O = 0.55 
in orbit around a G5 star. The linestyles correspond to log(g) = 2.3 (solid 
line), 3.0 (dashed line), 4.0 (dot-dashed line), 5.0 (dotted line). Left panel: 
PT-structures. Top right panel: PP-structures with pressure rescaled by 
^Q2.3-iog(g) figPf panel: Fraction of absorbed stellar flux as a function 

of pressure. Bottom right panel: Fraction of absorbed stellar flux as a function 
of rescaled pressure. 




Figure 10. Emission spectra as a function of surface gravity for a PefF = 1250 
K, [Fe/H] = 1 planet with C/O = 0.55 {upper panel) and C/O = 1.05 {lower 
panel) in orbit around a G5 star. The spectra are shown for log(g) = 2.3 (blue 
lines), log(g) = 3 (purple lines), log(g) = 4 (red lines) and log(g) = 5 (orange 
lines). The colored bars indicate the position of the absorption maxima of 
various species. The black line shows the blackbody flux at the atmosphere’s 
effective temperature. 


equilibrium one obtains the following simple relation between 
the optical depth r and the pressure P 

r=-P, (7) 

where k is the gray opacity and g is the gravitational accelera¬ 
tion (taken to be constant). In this case, changing the gravita¬ 
tional acceleration will conserve the temperature structure as a 
function of r, as r is the effective spatial coordinate for the ra¬ 
diation field. The mapping from r to P, however, will change, 
resulting in locations of a given optical depth and temperature 
to move to larger pressure values when g is increased. This 
is equivalent to saying that the location of the planetary atmo¬ 
spheric photosphere moves in terms of pressure if the surface 
gravity is changed. 

Thus, when plotting the PT -structures as a function of plan¬ 
etary gravitational acceleration, as can be seen in Figure 9, 
one notices that at higher log(g) the temperature structure 
appears to be shifted to larger pressures when comparing to 
cases with lower log(g). For demonstration purposes we show 
the RT-structures up to 10“^^ bar. Note, however, that we 
only calculate the opacities down to pressures of 10 “^ bar 
and adopt the 10 “^ bar values at all smaller pressures, i.e. 
k(P < 10“^ bar) = k(P = 10“^ bar). The RT-structures for 
pressures below 10 “^ bar are not necessarily unphysical, how¬ 
ever (see Section 4.3 for a discussion). The “highest altitude 
inversion” visible in this plot for pressures much smaller than 
10 “^ bar is due to the heating by the alkali line cores. 

In the top right panel of Figure 9 we show the PT -structures 
once more. In this case we have re-scaled the pressures in 
PT -structures with log(g) higher than 2.3 (which is the lowest 
log(g) value we consider) with To first order, his 

should counterbalance the pressure shift of the temperature 
structure induced by gravity when compared to the log(g) = 
2.3 case. However, as the opacities are non-gray and varying 
vertically we expect differences. Nonetheless, the resulting 
PT -structures lie on top of each other quite well. 

When comparing in greater detail one notices that the deep 
isothermal regions (at ~ 1-100 bars) are at higher tempera¬ 


tures for lower log(g). Here the pressure dependence of the 
opacity comes into play: for lower log(g) values the stellar 
light is absorbed at lower pressures, where the atomic and 
molecular lines are less broadened. This results in the stel¬ 
lar light being able to penetrate deeper in terms of rescaled 
pressure when comparing to high log(g) atmospheres. This 
means that more stellar light reaches regions of the atmo¬ 
sphere which are optically thick in the near-infrared, which 
does, in turn, heat up the atmosphere deep in these IR opti¬ 
cally thick regions. 

In the middle and bottom panel on the right side of Figure 9 
we show the fraction of the absorbed stellar flux with respect 
to the stellar flux at the top of the atmosphere. The middle 
panel shows this fraction as a function of pressure, the bottom 
panel shows this fraction as a function of rescaled pressure. 
One sees that the stellar light is able to penetrate deeper in 
terms of rescaled pressure in the case of low log(g). 

In Figure 9 we have shown an oxygen-dominated atmo¬ 
sphere, where the abundance of the main coolant and ab¬ 
sorber, H 2 O, is roughly independent of pressure. For carbon 
rich atmospheres the pressure dependent abundances of H 2 O, 
CH 4 and HCN might play a role in addition to the pressure 
shift of the temperature structures. 

In order to test that our above observations for the oxygen 
rich atmosphere are not caused by pressure and temperature 
dependent chemistry effects, we calculated self-consistent 
structures with vertically constant abundances of molecules 
and varied the surface gravity. We found the same behavior of 
the structures as described above, verifying that the pressure 
dependent line wing strengths are responsible. 

5.3.2. Spectra 

In Figure 10 we show the emission spectra of atmospheres 
with varying surface gravity for a planet with Reff = 1750 K, 
and [Fe/H] = 1 in orbit around a G5 host star, again for two 
different C/O ratios (0.55, 1.05). As mentioned above, a vari¬ 
ation in the surface gravity rescales the temperatures profiles 
in terms of pressure. We also found that the deep isothermal 
regions are hotter for the lower surface gravity cases, because 
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Figure 11. Atmospheric PT-structures and mass fractions as a function of 
metallicity for log(^) = 3 planets around a G5 star. The left panel shows 
the PP-structures of the cases with C/O = 0.55, Peff = 1500 K planets, the 
middle panel shows the cases with C/O = 0.95, Peff = 2250 K. The right 
panel shows the mass fractions of CH 4 (black lines), HCN (purple lines) and 
H 2 O (green lines) divided by the alkali mass fraction for the PP-structures 
shown in the middle panel. The different line styles in all panels stand for 
different metallicities: [Fe/H] = -0.5 (solid line), 0.0 (dashed line), 0.5 (dot- 
dashed line), 1.0 (double dotted dashed line), 2.0 (dotted line). 

the insolation can probe deeper into the atmosphere. In the 
pressure rescaled PT -structures (see upper right panel of Fig¬ 
ure 9) one can see that above the isothermal region the atmo¬ 
spheres of planets with higher surface gravity are hotter for a 
given rescaled pressure: The photosphere is located at higher 
pressures for a higher surface gravity. It is therefore less trans¬ 
parent, due to the line wing pressure broadening. In order 
to radiate away the required amount of flux the temperature 
therefore needs to be higher. The flux in the absorption fea¬ 
tures then originates in hotter regions, making the absorption 
troughs shallower in the C/O = 0.55 case. This behavior was 
verifled by the atmospheric structures with vertically constant 
molecular abundances as well. 

In the C/O = 1.05 case the same behavior can be seen, ex¬ 
cept for the atmospheres with the highest log(g), which shows 
emission features. Here the stellar light is absorbed over nar¬ 
rower and higher rescaled pressure ranges because the alkali 
line wings are much broader (the light is absorbed at higher 
actual pressure). The atmospheric cooling ability, however, is 
largely independent of pressure, because the emission of light 
depends on the Planck opacity /cp and dK^IdP = 0, if the pres¬ 
sure dependence of the chemistry is omitted.This causes the 
atmosphere at highest log(g) to develop an inversion. 

5.4. Metallicity dependence of the atmospheres 
5.4.1. PT-structures 

The influence of the metallicity on the PT-structures at 
low C/O ratios is as one would expect: An increased [Fe/H] 
value in atmospheres leads to higher temperatures in the deep 
isothermal part of the atmosphere in the cases where no in¬ 
versions are observed: the temperature structure is scaled to 
lower pressures as the metallicity increases, as a higher opti¬ 
cal depth is reached earlier in the atmosphere. The stellar light 
can penetrate deeper than suggested by a simple pressure scal¬ 
ing, however: the pressure dependent line wings are weaker 
(as the atmospheric structures shift to smaller pressures for 
higher metallicities). This increases the temperature of the at¬ 
mospheres in the deep isothermal regions at 1-100 bars (see 




A (/iin) 


Figure 12. Emission spectra as a function of metallicity for a = 1750 
K, [Fe/H] = 1 planet with C/O = 0.55 {upper panel) and C/O = 1.05 {lower 
panel) in orbit around a G5 star. The spectra are shown for [Fe/H] = -0.5 
(cyan lines), [Fe/H] = 0.0 (blue lines), [Fe/H] = 0.5 (purple lines), [Fe/H] 
= 1 (red lines) and [Fe/H] = 2 (orange lines). The colored bars indicate the 
position of the absorption maxima of various species. The black line shows 
the blackbody flux at the atmosphere’s effective temperature. 


left panel of Figure 11), just like it did for low surface grav¬ 
ities studied in Section 5.3. Similar to the test carried out 
for varying surface gravities in Section 5.3 we calculated test 
atmospheres with vertically constant molecular abundances, 
scaling the abundances by different factors for different struc¬ 
tures, mimicking variations in metallicity without having to 
deal with effects introduced by chemistry. These calculations 
showed the same behavior as the nominal calculations when 
varying the metallicity. 

In the case of FT-structures with C/O ~ 1, which have in¬ 
versions, the inversion temperature increases and the region 
directly beneath (i.e. at higher pressure) the inversions has 
a lower temperature if the metallicity increases (see middle 
panel of Figure 11). It is, at first, not evident why this should 
happen, because if all the metal atomic abundances scale with 
2 Q[Fe/H] would expect the same for the resulting molecular 
abundances and opacities, and therefore the heating vs. cool¬ 
ing ability of the atmosphere should stay the same. This inter¬ 
pretation is consistent with the analytical double-gray atmo¬ 
spheric models as published, e.g., by Guillot (2010); Hansen 
(2008); Thomas and Stamnes (2002), where the inversion 
temperature should stay constant unless the ratio 

r=^ (8) 

kiR 


changes, where /Cyis and a:ir are the mean opacities in the visual 
and IR wavelengths in the atmosphere. The behavior we see 
in the atmospheres suggests that 


dy 

d[Fe/H] 


> 0 , 


(9) 


which should only be possible if aTvIs and /cir (and the molecu¬ 
lar abundances giving rise to these opacitites) do not just sim¬ 
ply scale linearly with metallicity. In order to test this we 
checked the abundances of the major absorbers and emitters 
as a function of metallicity throughout the atmospheres for the 
PT -structures shown in the middle panel of Figure 11 . Indeed 
we found that the ratios of mass fractions 2fH2o/2f Alkali and 
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Figure 13. Atmospheric emission spectra of planets in orbit around a G5 star with = 1000 K. Left panel: Planets with log(g) = 4, [Fe/H] = -0.5. Right 
panel: Planets with log(g) = 2.3, [Fe/H] = 2. Upper subpanels: Emission spectra as a function of wavelength for planets with C/O = 0.55 (blue solid line) 
and C/O =1.12 (red solid line). The absorption bands of dominant absorbers are indicated by the colored bars below the spectra. Lower subpanels: Emission 
spectra as a function of wavelength (x-axis) and C/O ratio (y-axis). The flux values are indicated as a color map. The red-white dashed horizontal line in the right 
panel indicates the C/O value where the atmosphere switches from being rich in water to being methane-rich. The corresponding C/O value of this transition is 
indicated in the plots. The red and blue horizontal lines indicate the C/O values of the wavelength dependent spectra shown in the upper subpanels. 


^ch 4 /^Alkali decreased when the metallicity was increased 
(see right panel of Figure 11). ^hcn/^A lkali increases, at 
the relevant temperatures already being the dominant carbon 
opacity carrier. However, the increase in Xhcn/^A lkali is ap¬ 
parently not enough to act as a counterweight compensating 
the loss of infrared opacity due to the lower /^Alkali- This 

leads to less efficient cooling as [Fe/H] increases. This abun¬ 
dance change is likely caused by the pressure dependence 
of the chemistry, as higher metallicities shift the temperature 
structure to smaller pressures, where, for carbon-dominated 
atmospheres, CH 4 and H 2 O are less abundant, while the HCN 
abundance increases. 

5.4.2. Spectra 

Analogous to the log(g) case an increase in metallicity (and 
thus opacity) can be regarded as a similar pressure rescaling, 
as we found for a gray atmosphere with vertically constant 
opacity k that r = KigP. As k is in the numerator, atmospheric 
structures with increased metallicity should behave similarly 
to structures with J^creased surface gravity, featuring a higher 
temperature in their isothermal regions, but a lower tempera¬ 
ture (as a function of rescaled pressure) in the higher atmo¬ 
sphere: Because the photosphere will be located at smaller 
pressures (in actual, non-rescaled pressure) for an increased 
metallicity, the line wings will be less strong (less pressure 
broadening). The atmosphere is therefore more transparent 
and cools better. In order to radiate away the imposed flux, 
the temperature in this more transparent photosphere needs 
thus to be decreased. The minima in the spectrum, stem¬ 
ming from the opacity maxima, i.e. the line’s Gauss-cores, 
will originate from the same region in terms of rescaled pres¬ 
sure. As these pressures are now at a lower temperature, this 
leads to deeper absorption troughs in the spectra. This can be 
seen in the upper panel of Figure 12 and was confirmed by the 
vertically constant molecular abundance calculations as well, 
when rescaling the abundances as described above. In sum¬ 
mary, more pronounced absorption troughs can mean either 
a lower surface gravity or a higher metallicity (see figures 12 
and 10 ). 


In the C/O = 1.05 case we can again draw on our studies 
of the PT -structures: we saw that for atmospheres with inver¬ 
sions, due to the chemistry involved, the cooling ability of the 
atmospheres relative to the heating by the alkalis decreases if 
the metallicity is increased (see middle and right panel of Fig¬ 
ure 11). The spectra shown in the lower panel of Figure 12, 
although they do not exhibit inversions, are consistent with 
these observations, showing absorption spectra which become 
more and more isothermal as the metallicity is increased. 

5.5. Low temperature atmospheres (T^fi < 1250 K) 

At low enough temperatures (Teff < 1250 K) HCN does not 
yet play a significant role for the atmospheric spectra. Ad¬ 
ditionally the left pointing arrow of the chemical reaction in 
equation ( 6 ) can still be of importance, meaning that H 2 O and 
CH 4 are significant carriers of C and O atoms. It is important 
to note, however, that the chemical equilibrium abundances 
are not only temperature, but also pressure dependent: at low 
temperatures and high pressures CH 4 will be abundant also 
in an oxygen-rich atmosphere, while H 2 O will be abundant in 
carbon-rich atmospheres. At low temperatures and low pres¬ 
sures CO will become increasingly important, such that the 
oxygen-rich atmospheres do not contain a lot of methane and 
the carbon rich ones do not contain a lot of water. 

As seen in the above discussions, [Fe/H] and log(g) can 
strongly influence to which pressure levels the optical depth- 
dependent temperature structure will be be scaled, as for a 
gray atmosphere it would hold that r = KjgP. Therefore, low 
metallicity atmospheres (causing a small k) at high surface 
gravities cause the temperature structure to be scaled to high 
pressures. In Figure 13 we show emission spectra of planets 
with Teff = 1000 K in orbit around a G5 star. The spectra 
are shown for C/O = 0.55 and 1.12 in the upper subpanels. 
Furthermore we indicate the positions of absorption features 
of H2O, CO2, K, Na, CO, CH4 and PH3 in the plots. The left 
panel shows the emission spectra for planets with log(g) = 4, 
[Fe/H] = -0.5. This means that here the surface gravity is high 
and the metallicity is low, causing the temperature structures 
to be scaled to high pressures. The right panel shows planets 
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Tgff = 1250 K, condensation = 1250 K, No condensation 
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Figure 14. Emission spectra as a function of wavelength (x-axis) and C/O ratio (j-axis) of planets with Teff = 1250 K, log(g) = 3, [Fe/H] = 1 in orbit around a 
G5 star. The flux values are indicated as a color map. The red-white dashed horizontal lines indicate the C/O values where the atmospheres switch from being 
rich in water to being methane-rich. The corresponding C/O value of this transition is indicated in the plots. Left panel: Nominal chemical model (including 
condensation), right panel: Chemical model without condensation. 


with log(^) = 2.3, [Fe/H] = 2, i.e. with low surface gravities 
and high metallicities, leading to temperature structures to be 
scaled to low pressures. In the lower subpanels we show color 
maps of emission spectra as a function of wavelength (x-axis) 
and C/O ratio (j-axis). 

In the right upper subpanel, one sees that the two spectra 
are very different, showing either water or methane features 
for the atmospheres with C/O =0.5 or 1.12, respectively. As 
described above, this is expected, corresponding to a low pres¬ 
sure scaling of the temperature structure and due to the pres¬ 
sure dependence of the CO-CH 4 -H 2 O chemistry. In the lower 
right subpanel there is an overall shift from H 2 O to CH 4 dom¬ 
inated spectra at C/O -0.73. 

As expected, in the left panel there is only little differ¬ 
ence between the oxygen-rich and carbon-rich case. Fur¬ 
ther, the lower left subpanel does not show any transition be¬ 
tween a water- and methane-dominated atmosphere, as both 
molecules are present in the atmospheres at all C/O ratios. 
Once more, this is expected, as in this case, i.e. for low metal¬ 
licity and high log(g) the photosphere of the atmosphere is 
scaled to high pressures, where the chemistry dictates that 
CO is not the major carbon and oxygen carrier, but instead 
CH 4 and H 2 O dominate, at least at the low atmospheric tem¬ 
peratures considered here. Therefore, although the CH 4 /H 2 O 
number ratio may change as a function of C/O, this change is 
not sufficient to affect the spectrum significantly. 

Therefore, the spectral appearance of a planet is not only 
given by the C/O ratio and the effective temperature but also 
by a factor 

yS=[Fe/H]-log(g), (10) 

which is a measure for the optical depth - pressure mapping 
in the atmospheres and gives insight to which pressure levels 
a given atmospheric temperature profile T(t) is scaled. We 
found that transitions between water- and methane rich atmo¬ 
spheres occur at yS > -4 or -3.5 for = 1000 K. For Teff 
= 1250 K we found that fS > -5.0, indicating that a tran¬ 
sition between water and methane dominated spectra should 
always be expected at these temperatures. However, values of 

close to this threshold should always exhibit some methane 
or water features, even if the atmosphere is water or methane 
dominated, respectively. 


5.5.1. C/O dependence with and without condensation 

For atmospheres with effective temperatures < 1750 K, the 
spectrally active parts of the atmosphere have temperatures 
low enough for the condensation of MgSiOs (for the tempera¬ 
ture dependent saturation vapor pressure of MgSiOs see, e.g., 
Ackerman and Marley 2001). 

Condensation of O in MgSiOs does not have a too strong 
effect on the spectra in the sense that they are either water or 
methane dominated at high enough temperatures, i.e. Teff > 
1000 K. It does shift the C/O ratio-dependent transition be¬ 
tween the 2 cases, however, as we detail below. Note that we 
do not include cloud opacities yet, so the presence of MgSiOs 
grains will not affect the radiation field. 

For atmospheres with C/O values in the vicinity to, but less 
than 1, the condensation of MgSiOs decreases the amount of 
oxygen available to form CO and H 2 O considerably. In turn 
the H 2 O features in the spectra will weaken and CH 4 can form 
in noticeably higher abundances as C atoms are more avail¬ 
able due to the lower amount of CO being formed. 

This results in shifting the transition from H 2 O to 
CH 4 /HCN dominated spectra from C/O = 0.92, which we ob¬ 
tain for atmospheres with Teff > 1750 K, to C/O = 0.73 which 
we obtain for < 1750 K, as we described in the previous 
section. 

In order to test this condensation dependance further we 
carried out atmospheric calculations at = 1250 K, neglect¬ 
ing the effect of condensation. 

A comparison of the resulting emission spectra as a func¬ 
tion of C/O for both cases (Tgff = 1250 K, with and with¬ 
out considering condensation) can be seen in Figure 14. We 
calculated atmospheres with C/O ratios spaced equidistantly 
between 0.35 and 1.4 using 100 grid points for both cases. 
The difference in location for the shift from water to methane 
dominated spectra, moving from C/O = 0.73 (condensation) 
to C/O = 0.92 (no condensation), is very prominent in these 
plots. 

To further verify this finding we plot the mass fractions of 
H 2 O, CH 4 , CO and MgSiOs in Figure 15 for a planetary at¬ 
mosphere with Teff = 1250 K, C/O = 0.8, log(g) = 3, [Fe/H] 
= 1 in orbit around a G5 star. The C/O value is chosen such 
that the atmosphere is water dominated in the model neglect¬ 
ing condensation, but it is methane dominated in our nominal 
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Figure 15. Mass fractions of components in the atmosphere of a planet with 
C/O = 0.8, TefF = 1250 K, log(^) = 3, [Fe/H] = 1 in orbit around a G5 star. The 
solid lines show the mass fractions of H 2 O (blue), CH 4 (red), CO (magenta) 
and MgSiOs (black) for our nominal model, including condensation, while 
the dashed lines show the results for an atmosphere without condensation. 



Temperature (K) Mass fraction 

Figure 16. Left panel: PT -structure of the atmosphere of a planet with Teff 
= 1500 K, log(g) = 4, [Fe/H] = 1, C/O = 0.95 in orbit around a G5 star. Right 
panel: Mass fractions of CH 4 (red solid line), H 2 O (blue solid line), HCN 
(orange solid line), SiC (black solid line) and MgSiOs (purple solid line) as a 
function of pressure for the PT -structure shown in the left panel. 


atmospheric model, which includes condensation. 

One clearly sees that for high pressures, where the tempera¬ 
tures are too high for MgSiOs to condense, the abundances of 
H 2 O, CH 4 and CO for both models are nearly identical. The 
small differences are due to differences in the PT -structures 
found for the 2 chemical models. For pressures smaller than 
10“^ bar, however, MgSiOs starts to condense, noticeably de¬ 
creasing the CO and H 2 O abundances. CH 4 becomes much 
more abundant than H 2 O, which is in contrast to the behavior 
of the model without condensation, where H 2 O stays more 
abundant than CH 4 throughout the atmosphere. 

We therefore conclude that the transition from water- to 
methane-rich spectra may happen at C/O ratios considerably 
smaller than 1 if the planetary effective temperature is not too 
high. Especially for retrieval analyses of planetary spectra, 
which measure the local gas C/O ratio in the spectrally active 
regions of the atmosphere, the above findings are relevant. If 
condensation is expected to occur, a result such as “C/O <1”, 
due to the absence of methane features, could actually indi¬ 
cate an even lower total (gas -f condensates) C/O ratio<0.7. 
If a given atmosphere were enriched in Mg and Si one would 
expect this effect to be even stronger, shifting the transition 
between carbon and oxygen rich spectra to even lower C/O 
ratios. 

Finally, we want to remind the reader of the simplifica¬ 
tions of our chemistry model, which does neither include set¬ 
tling nor properly accounts for the effects of homogeneous 
or heterogeneous nucleation (see Section 4.2). Furthermore 
the absence of quenching in our models might be problematic 
if the timescales for condensation and chemistry in general 
are longer than the vertical eddy-diffusion timescales. Never¬ 
theless, similar results have been found with much more so¬ 
phisticated condensation models: Helling et al. (2014) were 
able to produce local C/O ~ 1-2 values in the gas phase for 
an atmosphere with a global C/O = 0.99 due to the conden¬ 
sation of O in dust species. Their model for condensation 
and cloud formation is much more complete and includes ho¬ 
mogeneous and heterogenous nucleation, settling, traces the 
growth and evaporation of grains, and considers more con¬ 
densable species than our model. 


Given these differences in condensation modeling it will 
be very important to reinvestigate our findings presented here 
with more sophisticated cloud models in the future. 

5.6. Intermediate temperature atmospheres (T^fi ~ 1500 K) 

At Teff = 1500 K the transition from oxygen to carbon dom¬ 
inated atmospheres is still at C/O = 0.73 as silicate conden¬ 
sation still takes place. Furthermore, the carbon dominated 
atmospheres show strong methane features, but HCN features 
start to emerge as well. 

5.6.1. Inversions and kinks at Teff = K and QO ~ 1 

At Teff = 1500 K and C/O ~ 1 condensation can lead 
to weak inversions and, occasionally, to kinks in the PT- 
structures. For C/O ratios ~ 1 the atmospheres are already 
carbon dominated. In general, the atmospheres are still too 
cool to contain enough HCN to efficiently radiate away the 
absorbed stellar energy, such that H 2 O and CH 4 are the main 
absorbers and the H 2 O-CH 4 -CO chemistry is important. At 
the intermediate atmospheric temperatures considered here, 
inversions are likely to occur because of the condensation of 
SiC. This results in a lower abundance of SiO, as less Si is 
available. The O atoms which are not bound in SiO anymore 
form more CO and thus decrease the C budget available to 
form CH 4 , therefore decreasing the atmosphere’s ability to 
cool. This effect can be further enhanced by the evaporation 
of MgSiOs in the inversion regions, which frees additional O 
to be put into CO, subsequently locking up more C atoms. 
As for the atmospheres which have inversions at C/O ~ 1 at 
higher effective temperatures, the inversions vanish for higher 
C/O ratios > 1: less oxygen is present to form CO in the first 
place. Therefore more CH 4 can be formed. 

At Teff = 1500 K kinks in the PT-structure can occur when 
condensation of MgSi 03 and the associated locking up of 
oxygen causes the cooling to become more strongly methane- 
dominated in a certain layer, whereas there is less methane 
present to cool in an adjacent, hotter layer in which there is 
less MgSiOs and more oxygen is available in the gas phase 
to form CO and water, locking up carbon and decreasing 
the methane abundance. The cooling ability of the methane- 
deprived layers is lower, leading to a strong temperature 
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Figure 17. Plots as shown in Figure 13, but for planets with [Fe/H] = 1, log(g) = 3 in orbit around a G5 star. Left panel: Teff = 1750 K, Right panel: Teff = 2250 

K. 


change from one layer to the next. Such kinks depend on 
the choice of the grid spacing and cell locations, such that 
they should not be treated as real physical phenomena but 
rather numerical artifacts. The corresponding structure files 
have been flagged with “_kink”. As inversions due to alkali 
heating and a low cooling ability are in general not seen in 
our results for M5 host stars, the kinks and inversions are not 
present for planets with M5 hosts. 

In Figure 16 we show an example for the kinks which are 
caused by the condensation. One clearly sees that the kinks in 
the PT -structures going to hotter temperatures are caused by 
the (partial) evaporation of MgSiOs, reducing the amounts of 
coolants such as CH 4 and HCN. 

5.7. High temperature atmospheres (T^fi > 1750 K) 

At high temperatures condensation processes do not play 
an important role anymore. Therefore, the transition between 
water and carbon-dominated spectra shifts from C/0 = 0.73 
to 0.92. Furthermore the carbon-rich atmospheres become 
more and more HCN dominated and CH 4 becomes less and 
less important as the temperature increases. As mentioned 
before, the chemistry is not only temperature but also pressure 
dependent, favoring HCN over CH 4 at high temperatures and 
low pressures. 


(or -2) in these atmospheres. For the larger (3 values the 
methane features fade away. 

Teff = 2250 K 

For Teff = 2250 K the atmospheres with C/0 > 1 are strongly 
HCN dominated. Only for low p values weak methane fea¬ 
tures are present. Furthermore more or less all atmospheres 
with C/0 ~ 1 have inversions if the spectral type of the host 
star is K or earlier. We show spectra of atmospheres with 
= 2250 K and varying C/0 ratios in the right panel of Figure 
17. 

Teff = 2500 K 

For Teff = 2500 K the atmospheres with C/0 > 1 are 
completely HCN dominated, and the methane features have 
vanished. All atmospheres with C/0 ~ 1 have inversions if 
the spectral type of the host star is K or earlier. 

6 . SUMMARY AND CONCLUSION 

In this work we present a systematic parameter study of hot 
Jupiter atmospheres. In addition to “classical” grid parameters 
such as metallicity, effective temperature and surface gravity 
we study the effects of the atmospheric C/0 ratio as well as 
the host star spectral type. We summarize the key findings of 
our study in Figure 18 and in the text below. 


Teif = 1750 K 

At Teff = 1750 K, we find that the higher the yS-factor 
(see Equation 10) of an atmosphere is, the more HCN 
dominates the spectrum. Methane features are visible for 
all ySs, however. Due to the chemistry, a low yS-factor allows 
for some presence of water in the carbon-rich atmospheres. 
Thus at low ySs we find a weak water absorption signa¬ 
ture imprinted on the rather opacity free region extending 
from 2.4-3 yum, which is bracketed by two CH 4 features. 
Because of the strong stellar irradiation the atmospheric 
structures at C/0 ~ 1 become either more isothermal or 
exhibit inversions. We show spectra of atmospheres with Teff 
= 1750 K and varying C/0 ratios in the left panel of Figure 17. 

Teff = 2000 K 

At even higher temperatures HCN becomes more dominant. 
Inversions at C/0 ~ 1 predominantly form for low p < -2.5 


• At low effective temperatures (Tgff < 1500 K) the 
atmospheres can be either water or methane dom¬ 
inated, but not always: \i p = [Fe/H] - log(g) is 
small, the spectra at < 1000 K are quite similar, 
showing both strong water and methane features. The 
optical depth (and hence temperature) versus pressure 
profile scales approximately with p. Hence, a given 
optical depth (temperature) is reached at high pressure 
when beta is low and vice versa. We want to remind 
the reader, however, that we neglect quenching, which 
could potentially alter the methane and water mixing 
ratios. For high pressures and low temperatures CH 4 
and H 2 O co-exist as the dominant oxygen and carbon 
opacity carriers, and are both visible in the spectrum. 
At p values above -4 to -3.5, the atmospheres are either 
water- or methane-dominated at Tgff = 1000 K. For at¬ 
mospheres with Teff = 1250 K the spectra look similar 
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only for the highest surface gravities (log(g) = 5) and 
lowest metallicities ([Fe/H] < 0), such that these atmo¬ 
spheres should be either water or methane dominated 
for most planets. 

At Teff < 1500 K the condensation of MgSiOs is a rel¬ 
evant effect at the local atmospheric temperatures. 

The condensation effectively lowers the amount of oxy¬ 
gen which can be put into CO and H 2 O, such that more 
carbon atoms are available to form CH 4 . As a result 
the atmospheres start to be methane dominated at 
C/O = 0.73. For higher temperatures MgSiOs can no 
longer condense, shifting the transition from oxygen to 
carbon dominated spectral signatures to C/O = 0.92. 

For planets with Teff > 1500 K and C/O ~ 1 host 
stars with spectral type earlier than M5 (we consider 
M5, K5, 05, F5) can lead to temperature inversions 
in the atmospheres. The reason for this is that under 
these circumstances all the main coolants of the atmo¬ 
sphere, H20, HCN, and CH 4 , are depleted, whereas the 
absorption of optical radiation by the alkali metals re¬ 
mains highly effective. For Tgff = 1500 K the condensa¬ 
tion of SiC can sufficiently lower the cooling ability for 
inversions to develop. At this effective temperature the 
condensation of MgSiOs can lead to kinks and numer¬ 
ical instabilities in the solutions for the RT-structure. 
For Teff = 2000 K all atmospheres with {3 < -2 io -2.5 
will exhibit inversions. For Teff ^ 2250 K all atmo¬ 
spheres with C/O ~ 1 exhibit inversions. 

The lower p = [Fe/H] - log(g )5 the more methane- 
dominated the spectra are at C/O ratios >1. At 

higher temperature and/or higher p values, such plan¬ 
ets have HCN-dominated spectra. In general we show 
the dominant absorbers as a function of temperature and 
C/O ratio in Figure 18. 

The host star spectral type is an important factor for 
the spectral appearance of the atmosphere. For plan¬ 
ets with C/O ~ 1 host stars of spectral type K or earlier 
can give rise to inversions if they are at small enough 
distances, whereas for M-type host stars inversions do 
not occur. Further, the later the host star spectral type, 
the more isothermal the planetary atmospheres become 
(if the C/O ratio is not ~ 1). This is because SED of the 
stellar irradiation becomes increasingly similar to the 
planetary radiation field. 

Planetary metallicity and surface gravity determine 
the location of the planetary photosphere. High sur¬ 
face gravities or low metallicities will shift it to larger 
pressures, whereas low surface gravities or high metal¬ 
licities shift it to low pressures. As the molecular and 
atomic line wing strength scales approximately linearly 
with pressure, for photospheres at low pressures the 
flux is originating at somewhat deeper layers (in terms 
of optical thickness), leading to somewhat cooler at¬ 
mospheric temperatures and hence deeper absorption 
troughs. The deep isothermal temperature increases in 
these cases, as the insolation can probe deeper into the 
atmospheres. Similar results for the surface gravity de¬ 
pendence of the absorption troughs have also been re¬ 
ported in Sudarsky et al. (2003). 
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Figure 18. Dominating IR absorbing/cooling species as a function of Teff and 
C/O. The red shaded region denotes carbon-dominated atmospheres, whereas 
the grey shaded region denotes oxygen-dominated atmospheres. The gray- 
hatched region denotes the temperature range where the atmospheric spectra 
can be dominated by CH 4 and H 2 O at the same time, independent of the C/O 
value. This occurs if [Fe/H] is low or log(g) is high. Within each region 
defined by the black solid lines the dominating IR absorbing/cooling species 
is indicated in the plot. The region in which inversions occur is shown in the 
plot as well. *Only host stars of type K and earlier can cause inversions. 
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It is interesting to see that at low temperatures the strength 
of methane or water features does not only depend on the 
C/O ratio, but also on the pressure level of the photosphere, 
which can be quantified using the the p factor. For higher 
temperatures the p factor plays a role as well, as it determines 
whether CH 4 or HCN dominates the spectra of carbon-rich 
atmospheres. Also the occurrence of an inversion at C/O ~ 1 
can be tied to the p factor, at least for the atmospheres with 
Teff ~ 2000 K. Therefore the p factor can be used as a third 
dimension to characterize the spectral appearance of an ex¬ 
oplanet, in addition to the effective temperature and the C/O 
ratio. 

Moreover, the fact that the transition from water to methane 
rich spectra shifts due to the condensation of silicates, which 
lock up oxygen, is important when carrying out retrieval 
analyses of planetary atmospheres. The C/O ratio is often 
measured by taking into account the abundances of only the 
gaseous carbon and oxygen carrying molecules. This can po¬ 
tentially overestimate the total (gas -f condensates) C/O ratio. 
It is important to note that our current condensation model is 
simplified, assuming instantaneous condensation once the sat¬ 
uration vapor pressure is exceeded and no settling or mixing 
of the cloud particles. It will therefore be very important to 
investigate this effect in the future in greater detail, using a 
more sophisticated condensation model. 

The fact that inversions can potentially occur at C/O ~ 1 
is interesting, as we did not require any additional absorbers 
such as TiO and VO, the absorption of stellar light by the al¬ 
kali atoms is sufficient. To further study the inversions it is 
necessary to obtain molecular line lists as complete as pos¬ 
sible as their occurrence is very strongly dependent on the 
atmospheric cooling ability. 

The grid of atmospheres presented in this work is made 
publicly available and can be found at the CDS. 


We thank the anonymous referee for a very thorough as¬ 
sessment of the manuscript and many comments that have 
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APPENDIX 

A FAST METHOD TO CALCULATE OPACITIES FROM LINE LISTS 

If one wants to calculate the line opacities of a given molecule on a given grid of wave number points one must, in principle, and 
if no line truncation is applied, calculate the line profile of every line at every wave number grid point. Even if one precomputes 
the opacities and tabulates them for later use it can still take a long time to calculate the total molecular cross-section as the 
calculations need to be carried out at high resolution. Using a fiducial resolution ofR= 10^ and considering a species with of the 
order of 10^ lines one easily ends up with ~ 10^^ line profile evaluations at just one given pressure and temperature. 

To speed up the calculations of line opacities the method explained below was developed and used for our opacity database 
calculations. In essence we calculate the line cores of every line at high resolution, while calculating the line wings far from the 
line core on a much coarser grid. Once the contribution of all the lines to the coarse grid has been calculated it is interpolated 
back to the fine grid. In detail we proceed as follows: 

Divide the total wave number grid into subgrids of 10,000 grid points. Then start to go through all these subgrids, which will 
be indexed by m: 

• Calculate all line opacities of lines that are lying within the subgrid m at all of its 10,000 wave number grid points. 

• Then iterate over all other subgrids (which are indexed by n). For the external lines in a given external subgrid with n ^ m 
do the following: 


1. If 7g > Jl (i-e. Gaussian width larger than Lorentz width): If the distance of the line to the subgrid border of m is 
smaller than /g7g. where /g is a factor that needs to be specified, then this line gets calculated at all of m’s 10,000 
subgrid points. Otherwise go to step 2. 

2. If 7l > 7g or the distance to the subgrid border of m is larger than /g7g- 

Consider the Lorentzprofile yliy^ + {x- aq)^). Far away from the line center its functional form is roughly yl{x- aq)^. 
The relative deviation a from this form is 


+ (x - xpf r y _ y 

y (x - xo)2 y"^ + {x- xo)^ 


(x - Xo)2 


(Al) 


I.e., if we want a maximum deviation of less than a from the above form, then we need that 


|a-A ol > —. 


(A2) 


If this is not fulfilled, then the line opacity is just calculated at all of m’s 10,000 subgrid points. It it is fulfilled go to 
step 3. 

3. For all lines within a given external subgrid n that fulfill the above inequation: Calculate the line strengths on a coarse 
subgrid of 10 points within the original subgrid m and add the results for all these lines up, then interpolate back to 
the 10,000 original grid points in m, using a powerlaw interpolation and a coordinate transformation (a simple shift). 

4. Move on to the next external subgrid ^ -r 1, go back to step 1. 


• Move on to the next subgrid m. 


The reason to interpolate back to the 10,000 subgrid points of m for all external subgrids n seperatly is the following. 
For a single line, far away from its line center, the line shape is roughly 


(/){x) = 


7 

(x - Xo)2 ‘ 


(A3) 


Thus for all transitions of strength S k„ which are contained in a subgrid n the total continuum line strength (i.e. the wing 
strength of a line far away from its line core) will be 


^ tot,n 




kn 


S k^ykn 
(a - XkJ^ ■ 


(A4) 


Seen from subgrid m all lines kn within a given subgrid n have roughly the same line center position (namely within subgrid n), 
thus one can do the coordinate transformation y„ = a - a„, where a„ is the position of subgrid n (e.g. the wave number at its 
center). This would yield 



S k^ykn 

(yn-^Xn- XkJ^ 


1 

yl 


k„ 


(A5) 
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Figure 19. Schematic drawing of the algorithm used to calculate the opacities. For the lines within subgrid m the opacities get calculated on the fine fiducial 
grid. If a subgrid n mis sufficiently far away (in this example n=m+2) the opacity of the lines in n are calculated on a coarse grid in m. The summed opacity 
values of the lines in grid n are then interpolated on the fine grid in m using a coordinate shifted powerlaw interpolation. 


Thus for every subgrid n one can do a coordinate transformation to y„ and finds that the coarse 10-point subgrid continuum of 
subgrid n seen in subgrid m should roughly behave like a powerlaw function with a powerlaw slope of ~ -2. This explains 
why a powerlaw interpolation in the coordinate y„ is the best thing to do when interpolating the coarse continuum of the lines in 
subgrid n back to the fine grid in subgrid m. It should also be stressed that it is better to do an interpolation that finds the effective 
powerlaw slope, rather than taking it to be -2 and using there will be slight deviations from this -2 powerlaw 

shape, as one knows that the line centers in subgrid n are close to Xn, but not exactly at An interpolation will mitigate this 
problem by finding a slightly different powerlaw shape and an overall coefficient for the function that slightly deviates from 
Y^kn^knJkn' ^^r contmua produced by lines that are outside of the total grid we use a linear interpolation, as we don’t actually 
check where these lines are sitting. In Figure 19 one can see a schematic drawing explaining the acceleration method introduced 
in this section. 

All our opacity calculations were carried out using the accelerated method on a grid with a point spacing of T/A/l = 10^. 
Additionally we performed a “classic” calculation on a reduced grid constructed by using every 1,000th fiducial wave number 
grid point. On this reduced grid we did not use the aforementioned acceleration method but calculated every line contribution at 
every wave number point. The high resolution result of the accelerated method was only kept if the maximum relative deviation 
at the points coinciding with the 1,000 times coarser test grid was smaller than 1 %. If it was bigger, fc was increased and a was 
decreased and the calculation was repeated. 

CORRELATED-K: GOING FROM O (c^) TO O (A) 

In the following we describe our method of combining the opacity tables of multiple species. Furthermore its implementation 
at different grid resolutions is explained. For a general review of the correlated-k method see, e.g., Marley and Robinson (2014). 


The ‘'classical” 0[c^^ case 


The commonly utilized method to combine the k-tables of multiple species is numerically quite expensive, as it is of order 
O where Ng is the number of grid points used in g-space (g is the cumulative opacity distribution function, see below) and 

Asp is the number of species. In this traditional method, the computation of the total opacity ATtot works as follows: In a spectral 
region of the frequency interval [v, v -l- Av] the transmission of light T through a layer of thickness AP which contains 2 spectrally 
active species is 


T = 


I 


v+Av 


exp 


Xi/(i(y) + X 2 /( 2 (y) . o 

- AP 

a 


A^ ’ 


(Bl) 
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where Xt and Kt are the mass fractions and opacities of the two species, a is the gravitational acceleration in the atmosphere and 
AP is the atmospheric layer thickness in units of pressure. For simplicity it is assumed that Xi and Kt are constant within the 
atmospheric layer. If one assumes the opacities of species 1 and 2 to be uncorrelated, i.e. 


ftot(Kl,K2) = Mki) ■ f2iK2) , 


(B2) 


where / are the opacity distribution functions, one can rewrite the transmission T as 


T = 


I 


■y+Ay j / 

^-XiKi(y)AP/a^ 

Ay 


f 


>y+Ay j / 

^-X2K2(v)AP/a 

Ay 


(B3) 


An opacity distribution function within a frequency interval [v, y+Av] is defined by f(K)dK being the fraction of the opacity values 
within [y, y + Ay] which lie between k and k + dA:. Going from frequency space to g-space, where g is the cumulative opacity 
distribution function (dg = f(K)dK), and approximating the integrals with sums yields 


Ns N, 

J Jexp 

j=l 


+ X 2 K 2 J 


AP 


^gAgj 


(B4) 


The combined total k-table of species 1 and 2 therefore has the opacity values 


which have to be weighted with 


^tot ,/7 ~ + X2K2J 

Agij = AgiAgj. 


(B5) 

(B6) 


As is commonly pointed out the number of operations that need to be carried out in order to combine the k-tables of multiple 

species is thus 0(Ng'^), which can make the consideration of multiple species computationally expensive (see, e.g., Marley and 
Robinson 2014; Lads and Oinas 1991). 


The O (N) case 

Algorithm used at a bin size ofA/AA = 1000 (RIOOO method) 

In order to combine the individual k-tables of all species of interest for finding the total k-table of an atmospheric layer we use 
a method that is computationally less expensive. Similarly to the “classical” approach, the method makes use of the assumption 
that the opacities are not correlated. The main idea is to iteratively combine the opacities of two species: The opacity of a real 
species and the effective opacity of a “help”-species. If the opacities of all species are uncorrelated, then the combined opacity 
of two species is not correlated with the opacity of any other remaining species. Furthermore the combined opacity of the two 
combined species can be treated as belonging to a new single species, which is the “help”-species. 

We thus proceed in the following way: For every species, within every Ay bin, we save the opacity distribution /c(g) on a grid 
of 30 points. The 30-point grid consists of two 15-point Gaussian grids ranging from 0 to 0.9 and from 0.9 to 1, respectively.^^ 
Now, when starting to construct the total opacity, the first two species 1 and 2 get combined according to equations (B5) and 
(B6). This results in 30 x 30 = 900 new values K\+ 2 ,ij which need to be sorted by size. Using the cumulative sum of the associated 
weights AgiAgj, where Agi and Agj are the respective Gauss-grid weights, we interpolate the result back to the original 30-point 
Gauss-grid. This newly obtained opacity K 1+2 is then iteratively combined with the remaining species’ opacities and results in the 
final opacity distribution /CtotCg). In a procedural notation the method can thus be described as 

Total opacity = X_1 * kappa(g) of species 1 

For all remaining species (i = 2 to N_sp) { 

Total opacity = combine(Total opacity, 

X_i*kappa(g) of species i) 

re-bin Total opacity to nominal g-grid 


} 

For notational convenience the method for combining the opacities as introduced in this section will be called “RIOOO method” 
in the following sections. The number of points used for combining two species’ opacities will be called A^rics- As explained in 
this section, the nominal value of Miies when working at a resolution of A/AA = 1000 is Arics = 30. 

This is not the same grid on which the radiative transport will be carried out on. The radiative transport grid consists of 20 points. The 30-point Gaussian 
grid is only used for the combination of the k-tables. 
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Algorithm used at a bin size of XjAX - 10 and 50 

In order to correctly describe the opacity distributions at small XjAX many g-grid points would need to be used, as especially 
at low pressures the opacity tables K(g) tend to be very sharply peaked at g values very close to 1. Therefore the RIOOO method 
would become numerically inefficient and cannot be used. 

However, once more the idea is to combine two species iteratively in order to obtain the total opacity. Moreover, we again 
make use of the assumption that the opacities are not correlated. In the X/AX =10 and 50 cases the spectral bins are 100 or 20 
times larger than in the X/AX = 1000 case. They therefore include many lines, and the assumption of uncorrelatedness should be 
valid to an even higher degree than in the X/AX = 1000 case. 

The idea to obtain the total opacity is the following: In principle the combination of two species could be accomplished by 
randomly sampling the 2 individual opacity distributions and taking the sum of the sampled values as a set of the combined 
opacity. In a numerically simplified version one could discretize the opacity distributions by providing a pre-sampled set of N 
opacity values and their corresponding weights Ag. 

The random sampling could then be approximated by randomly drawing values from the opacity sets of each species and 
adding them, taking into account their weights at the same time. If one would sample continuous values from a distribution, it is 
possible to sample values from within a given interval multiple times. Thus, if a discretized opacity value has been drawn from 
the opacity set it must in principle not be excluded from being drawn in any of the next sampling steps. 

The discretization is carried out in the following way in our method: For every species we divide the K(g) table of every species 
into two sets. The first set contains K(g) values with g < gbord- The g-coordinates are located at the centers of grid cells defined 
by + 1 grid borders spaced equidistantly between g = 0 and g = gbord- The second set contains K(g) values with g > gbord- 
These g values are located at the centers of Np grid cells defined by + 1 grid borders spaced equidistantly between g = gbord 
and g = 1. We chose gbord = 0.985 and Np = 128. The Np values of a species with g < gbord will in the following be called a:iow 
and the Np values with g > gbord will be called /Chigh- Additionally, for every species, we save the lowest and highest opacity value 
within the frequency bin, corresponding to the g = 0 and g = 1 opacity values, /ciow describes the low g, continuum properties of 
the species’ opacity, while AThigh describes the high g, line core properties of the species’ opacity. 

Returning to sampling values from 2 species, the probability of sampling and combining 2 values stemming from the respective 
g < gbord-i'egions is The probability for combining 2 values from the g < gbord-i'egion of species 1 and the g > gbord-i'egion 
of species 2 is gbord • (1 - ^bord) etc... 

To speed up sampling, we now assume that once an opacity value of a given species has been drawn, it cannot be drawn again 
(we will return to the validity of this approach below). 

In order to approximate the sampling process of the combined opacity distribution function of two species, we then construct 
a 4Np X 2 matrix K containing the various possible combinations of a:iow and AThigh of both species, weighted by how common 
these combinations would be in a random sampling process of both species’ opacities. When sampling points from species 1 and 
combining them with sampled points from species 2 the assumption that a given value can not be redrawn allows for a simple 
shuffling in the sampling process: 


' Xi • shuffie(/^ijow) + ^2 • shuffie(/^2,iow) 
Xi • shuffie(/^ijow) + X2 • shuffie(/^2,high) 
Xi • shuffie(A:i,high) + ^2 • shuffie(A:2,iow) 
• shuffle(/^i,high) + ^2 • shuffie(/C2,high) 


sLd \ 

Np 

^bord(l ^bord) 

Np 

gbord (1 gbord) 

(1 gbord) 

Np / 


(B7) 


The first column of K represents the sampled values of the new combined opacity, the second column gives the weight of each 
sampled value, similar to the AgiAg 2 weights in the classical method described in Section B.l. We then sort the lines of the 
matrix K by the values in the first column. After this we construct a vector y of length 4Np with the entries (starting at m = 2) 


ym ~ ym-l X' 


k(m-l),2 + km,2 
2 


(B 8 ) 


and yi = ki^/^. The second column of K is then replaced with y. After this, the first column of K contains the newly sampled 
A:tot(g) values of the combined opacity of species 1 and 2, the second column contains the corresponding g values. Using A:tot(0) = 
Xia:i( 0 ) + X2K2(S4) and A:tot(l) = XiK\{\) + X 2 a: 2 ( 1 ) the total opacity can then be interpolated to the Np low-g and Np high-g values 
to yield the final result. The resulting opacity is then ready for being combined with the opacity of the next species. In order to 
shuffle the opacities we use the Knuth-shuffle algorithm, which is of order 0{Np). 

The assumption of not being able to draw a given opacity value more than once is obviously not correct. However, it has been 
found to not affect the quality of our results. From the above we see that in every combination step one needs to sort 4Np values. 
In the RIOOO method we would have the same computational costs when storing A^Ries = 2 ^JNp opacity points per species, 
losing resolution when comparing to the 2Np points we use in the method introduced here. Furthermore, the results of the RIOOO 
method, at the same computational cost, turn out to be much worse, both when comparing to the actual shape of the wanted total 
opacity distribution as well as when comparing 


1 

A^ 


pv+Nv 

1 ' 


ATy/dv' 




(B9) 
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— New method - - = 30 = 256 

Figure 20. Comparison of the different methods to combine the K(g) tables of different species. Upper panel: Opacity of water (red solid line), methane (red 
dashed line) and ammonia (red dotted line) as a function of g. The total K{g) obtained from adding the opacities in frequency space is shown as a red thick solid 
line. The results when using the RIOOO method with Nrics = 30 points and Nrics = 256 points are shown as black dashed and dotted lines, respectively. The 
result when using the new method introduced in this section is shown as a black solid line. Lower panel: Relative error of the three methods compared to the 
correct solution: Nrics = 30 points (dashed line), Nrics = 256 points (dotted line), new method (thick solid line). 

for both methods.The error of our method is in the range of %, whereas the error of the RIOOO method at the same computa¬ 
tional cost is in the range of tens of %. Comparing the results of the new method with results of the RIOOO method when taking 
A^RieS = i-G- the same number of points in both cases, yields slightly better results for the RIOOO method. However the 
numerical costs for the RIOOO method are 0{AN^), while they are 0{ANp) in the new method presented here. The reason for the 
RIOOO method at the fiducial resolution A^Ries = 30 to fail here is that we consider 20-100 more points per wavelength bin. This 
requires a higher resolution when trying to resolve the actual opacity distribution function. 

In Figure 20 one can see an example calculation from combining the opacities of water, methane and ammonia in the wave¬ 
length range going from 6.64 to 7.34 yum. The K{g) distributions of the individual species are shown in the plot. All species are 
contributing approximately equally strong to the total opacity in this example and have lines in the wavelength region of inter¬ 
est. Therefore this case represents something like a worst-case scenario, as our method is the most accurate when one species 
dominates or the other species only contribute via a their line continua. We plot the correct total /c(g) distribution, obtained when 
adding the opacities in frequency space first, as well as the results obtained from using the RIOOO method and the result from 
using the new method introduced in this section. The g-grid used for the RIOOO method was chosen to have g values following 
a distribution cx dlogA:/dg in order to trace strong changes in the opacity distributions. One sees that our new method is never 
worse in accuracy than the RIOOO method which even has a little higher computational cost (A^rics = 30), and usually has an 
relative error which is an order of magnitude smaller. The error of the A^Ries = 256 results is an order of magnitude smaller than 
the Arics = 30 result. 

Finally we note that our spectral calculations using the above efficient method at T/A/l = 10 do not deviate by more than 5 % 
(and usually less) in wavelength regions of appreciable flux when comparing to the rebinned 4/AT =10^ line-by-line calculations 
(see Section 2.4.1, Figure 2). This is a deviation commonly stated for correlated-k (see, e.g., Fu and Liou 1992; Lacis and Oinas 
1991). The strength of the new method reported here is to be numerically efficient, while conserving the opacity information at a 
high level of detail. 


USE OF THE VARIABLE EDDINGTON FACTOR METHOD TO FIND THE TEMPERATURE 

Basic equations 

For the moment equation based approach of solving for the PT -structure we first need the equation of radiative transport, 
neglecting scattering processes for now: 

n • V/y(x, n) = -ay{x) [/y(x, n) - 5'y(x, n)]. (Cl) 

In our case, the source function Sy is simply the Planck function, 

Sy(x,n) = By(n (C2) 


where T = T(x). 

We now make the plane-parallel assumption, which states that any spatially varying quantity can only vary in the vertical direction 
z. We chose z to increase towards the upper layers of the atmosphere. The equation of radiative transport then transforms to 

n^ly{z,n,^) = -ay(z) [Iy(z,Sy(Z,IJ., (/>)], (C3) 

az 

We will need to evaluate Eq. (B9) when computing the Planck mean opacity in the temperature calculation. 
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where fi = cos(^) and 6 being the angle between the vertical and the direction of the ray. 0 is the polar angle around the z-axis. 
Note that S is independent of both yu and 0 when it is equal to the Planck function. 

The zeroth, first and second radiative moments are defined as 


[/yCx), Hv(x), ^v(x)] = ^ ^ ^v(x, n) [1, n, nn] dfi. 


(C4) 


In plane-parallel geometry and rotational symmetry around the z-axis (i.e. no ^-dependence), only the z-component of the first 
moment H and only the zz-component of second moment K are unequal to 0 and one can define 


ff(z) = H,(z), 

K(z) = K,,{z\ 

where the v subscript has been omitted. The definition of the three plane-parallel moments then is 


[/, 


1 

■{z),Hy{z), Ky{z)] = 2 J 


For radiation emanating from a small solid angle (while keeping the z-only spatial dependancy) one finds that 

[Jy{z\Hy{z),Ky{z)] = /,,y (Z, jU, ) [ 1 , jU. , /T« ] , 


(C5) 

(C6) 


(C7) 


(C8) 


where = cos(^=,) and 6^ being the angle between the vertical vector and the vector pointing in direction Hy{z) and Ky{z) 
are, once more, the z- and zz-component of H and K. If the radiation emanates from a star of radius at distance J, where 
d then 




(C9) 


Integration of Eq. (C3) over the whole solid angle yields 


^Hy{z) = -ay(z) [Jy(z) - By{z)] , (CIO) 

dz 

where we used that the source function is supposed to be the Planck function. Note that this equation holds independently of the 
fact whether there is a 0-dependance in the radiation field or not as long as the definition H{z) = H^i^z) is used. It can thus also be 
used for the radiation emanating from a small solid angle. Multiplying Eq. (C3) by yu and integrating over the whole solid angle 
again yields 

^Ky{z) = -ay{z)Hy{z), (C11) 

dz 

where the isotropy of the Planck function was used. Equations (C7), (C8), (CIO) and (Cll) are the equations of interest for the 
task of finding the PT structure. 


Solution of the PT-structure problem 

The method explained below is based on the method used for protoplanetary disks introduced in Dullemond et al. (2002). 

A useful spatial coordinate for the PT-structure calculation is the pressure P, rather than the height z: At the top of the 
atmosphere we have z ^ oo and the starting point of z = 0 can be chosen arbitrarily. In contrast to the pressure, where we have 
a well defined value at the top of the atmosphere, namely P = 0. Furthermore the use of the pressure instead of some arbitrary 
height z in the atmosphere makes the density drop out of all equations of interest. To eliminate the vertical height z from the 
equations, we use the equation of hydrostatic equilibrium 


dP = -pgdz, (C12) 

where p is the density and g is the gravitational acceleration, which is taken to be constant throughout the atmosphere. Further¬ 
more we will use that 

Qfy(z) = pKy{z), (Cl 3) 

where Ky is the monochromatic opacity. 

The optical depth Ty then relates to the height z as 

dry = -pKy&z. (C14) 


dry = — dP. 
g 


This yields 


(C15) 
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Incident stellar irradiation 

For simplicity we assume the stellar irradiation field to be a black body in this section. In the implemented version of our code 
we are able to use either a blackbody or an actual stellar spectrum. In the latter case gets replaced with the stellar spectrum 
appropriate for a main sequence star at a given effective temperature. 

First we start with the stellar light shining at the atmosphere of the planet. The stellar effective temperature shall be T^. If one 
then defines the irradiation temperature as 

(R 

Tirr = (^) n (C16) 

then equations (C8) and (C9), together with Iy^^{P = 0) = By{T^) yield for the frequency integrated first moment in z- direction 
that 

lu.a-r' 

H,{P = 0) = -- 


where it was used that 


f 


irr 

471 

(C17) 

^ rj.4 

^ * • 

71 

(Cl 8) 


The negative sign implies that the radiation enters the planet, rather than leaving it. Furthermore we make the so-called two 
stream approximation, assuming that the stellar irradiation is in the optical wavelengths, while the radiation field inside the planet 
is in the IR-wavelengths due to the lower temperature of the planetary atmosphere. Any emission processes in the atmosphere at 
the stellar irradiation wavelengths are thus neglected and the corresponding source term in the equations of interest are neglected, 
leading to 

= jy,* (C19) 

dTv 


and 


dTv 


-Ky 


Hy 


(C20) 


where equations (CIO), (Cll) and (C14) were used. Using Eq. (C8) to see that Ky^, = pzJy^* yields 


d^ 1 


PZ 


For attenuation in the atmosphere one then finds that 


H,(P) 


0 


Hy^^iP = 0)e 


-Ty/fl. 


dv 


with 


1 

= - Ky(P')dP'. 

g Jo 


(C21) 


(C22) 


(C23) 


The important equations from this section are equations (C22) and (C23). 

In the case of taking the day side or global average of the stellar radiation we assume the stellar irradiation to be isotropic. In 
this case we use that 

4(y, P = 0) = -4i7.(y, P = 0), (C24) 


which follows from 


■i£ 


pP(y, P = 0)dp 
phiy, P = 0)dp , 


(C25) 


where we used that the stellar light only shines downward in the last line. Further assuming that at the top of the atmosphere h is 
independent of p (isotropy) leads to the desired result. We then carry out a full angle and frequency dependent radiative transport 
calculation for the stellar intensity 4, assuming only attenuation. From this we can calculate the stellar flux in every layer. 


Planetary radiation field 

The total net flux leaving the planet is supposed to be crT^^^. As the planet receives -p^crT^^, the wavelength integrated flux 
coming from within the planet at IR wavelengths must be 


471 


4n 


H(P = 0) = 


(C26) 
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The total flux is 




(C27) 


i.e. 


HtotiP = 0) 


471 


(C28) 


as required. As there are no sinks or sources of energy for the radiation held in the steady state equilibrium case we know that 


dP 


0 . 


Together with equations (C22), (C27) and (C28) this yields 

o-Tf 


crT^ 

H(P) = - Hy,(P = 

47r Jo 


(C29) 


(C30) 


For the solution one uses the wavelength dependent opacities of the previous full radiative transfer step to calculate the attenuation 
of the stellar light. 

As the next step we deflne the /y-averaged Eddington factor f as 


/ 


_ 1 r 

~~J Jo 


fyJydV 


K 

7 ’ 


(C31) 


where K and J are the wavelength integrated moments of zeroth and first order and /y = KyjJy. Eq. (Cll) then yields, together 
with / and equations (C12) and (C13): 

d 1 r 

dR g Jo 


KyHydV 


(C32) 


and Anally 


(Z*^) - -krH, 
dP g 


where kr is the Hy averaged opacity. For the solution of J one takes 

J{P = 0) = -H(P = 0), 

and uses /(R) and kr{P) of the previous full RT step and the results of Eq. (C30) to integrate Eq. (C33) from R = 0 to the 
pressure of interest. 


(C33) 


(C34) 


Finding the temperature 

Once one has obtained /(R) of the planetary radiation held one can use the wavelength integrated version of Eq. (CIO), noting 
that E^tot = constant vertically. This yields 


-T\p(T)-KjJ- 

71 


J ^oo 

0 


JyAV = 0, 


(C35) 


where a:/ is coming from the previous full radiative transfer step. As one finds from a similar analysis as performed for Hy ,, that 

4JP) = (C36) 

and as Eq. (Cl9) gives that 

’ (C37) 


one Anally gets 


4,(P = 0) =- HyJP = ^) 


^OO 

-RV(R) -kjJf Ky—HyJP = = 0, 

TT Jo R* ’ 


(C38) 


which has to be solved for T to And the temperature at pressure R for the next iteration step. In the code this was done by applying 
a zbrent root flnding algorithm taken from numerical recipes (Press et al. 1992). Furthermore, every 20th iteration step we evolve 
the temperature structure by using a Ng-accelerationn (Ng 1974) applied on R^. 

From the corr-k full radiative transfer step we thus need the opacities of the previous iteration step for calculating the attenuation 
of the stellar light, the /y-averaged Eddington factor, if/, kr and a:/ as well as Kp to And the pressure temperature structure. 
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Treatment of convection 


After the radiative structure of the atmosphere has converged we switch on convection in our code. 

During the moment solution of the temperature, the radiative temperature profile is solved from top to bottom (starting at low 
Pq, typically Pq = 10“^^ bar). We check in each layer i whether it should be convective or not by comparing the effective radiative 
temperature gradient 


^rad — 


Tj - Tj.i 

Pi-Pi-1 


l Pi^Pi-i \ 


(C39) 


with Vad = (r 2 - l)/r 2 , where 


r2 = 



p 


Xt 


cppT Xp 


(C40) 


with P being the pressure, T the temperature, p the density, cp the specific heat capacity, Xt = {d\ogPld\ogT)p and Xp = 
{dlogP!d\ogp)T (see, e.g., Hansen et al. 2004). All required quantities can be obtained from the equilibrium chemistry code CEA. 
We evaluate r 2 as r 2 = (r 2 ,/ + r2,i-i)/2 on our grid. 

We employ the Schwarzschild criterion, such that if Vj-ad > Vad, we adjust the temperature in layer i to be 


^ Pi.i^Pi(2r2-l) 

‘ '■‘>, + p,_i(2r2-i) ■ 


(C41) 


As the energy in a convective layer is not transported by radiation anymore, the integration of J via Eq. (C33) is not possible 
in this layer. However, in order to be able to discriminate between radiative and convective energy transport in layers lying below 
a current convective layer we need to compare to the radiative temperature in deeper layers. For this we need to continue to 
computation of J down to deeper layers. 

We thus chose the approach that in a convective layer i during the ^-th iteration = withofpi = 

with being the mean intensity taken from the full angle and frequency dependent radiative transfer step of the previous 
iteration. The superscripts indicate the iteration number from which the respective quantity is used. We chose this approach as in 
the case of very efficient convection (Viayer ^ Vad) the atmospheric layer should be optically thick, i.e. // ^ B(Ti). When going 
to the next layer / + 1 we radiatively integrate J to this next layer using Eq. (C33) and compare the resulting Vj-ad with Vad again. 

As in Marley et al. (1996); Burrows et al. (1997) we only allow a limited number atmospheric layers to be changed 
to convective energy transport every iteration. This is done to allow the atmospheric structure to adapt to the introduc¬ 
tion of convective layers. In Marley et al. (1996); Burrows et al. (1997) only 1 layer per iteration is allowed to change. 
We allowed for the change of 2 layers per iteration, because sometimes a layer on the brink to being convectively unstable 
will switch back and forth between being radiative or convective, preventing the overall convergence of the atmospheric structure. 
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